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ABSTRACT 


New  micromechanical  models  for  the  prediction  of  particulate  composite  mechan¬ 
ical  behavior  have  been  developed.  The  models  use  an  energy  balance  concept  to  account 
for  nonlinear  behavior  due  to  particle  debonding  and  incorporate  a  composite  modulus 
prediction  routine  based  on  an  improved  Mori- Tanaka  method.  This  method  permits  par¬ 
ticle  interaction  effects  to  be  taken  into  account  and  allows  the  stiffness  matrix  for  voids 
or  vacuoles  to  be  explicitly  stated  in  the  model.  To  demonstrate  the  characteristics  of 
the  improved  Mori-Tanaka  method,  comparisons  with  2-phase  and  3-phase  modulus  data 
were  made.  The  micromechanical  models  developed  for  void  and  vacuole  formation  were 
evaluated  against  available  mechanical  behavior  data.  Comparisons  showed  that  the  model 
derived  for  vacuole  formation  predicted  the  mechanical  behavior  correctly  for  two  types  of 
model  composites  which  contained  inclusion  volume  fractions  ranging  from  0.2  to  0.5.  The 
inability  of  the  models  to  predict  the  initial  stress-strain  behavior  of  a  composite  containing 
a  volume  fraction  of  0.22  of  well-bonded  particles  suggests  that  the  nonlinear  matrix  effects 
need  to  be  included  in  the  present  model  formulation. 


RESUME 


Les  nouveaux  modeles  micromecaniques  pour  la  prediction  de  comportement  des 
composites  charges  ont  ete  developpes.  Les  modeles  utilisent  un  concept  d’energie  pour  tenir 
compte  de  la  perte  due  aux  decollements  de  particules.  Les  modeles  emploient  aussi  une 
technique  crees  par  Mori  et  Tanaka  pour  calculer  le  module  de  composite.  Grace  a  cette 
methode,  on  peut  inclure  dans  nos  equations  les  effets  d’interaction  entre  les  particules 
et  les  effets  de  types  de  decollement  soit  les  vides  ou  les  vacuoles.  Pour  demontrer  les 
caracteristiques  de  cette  nouvelle  methode  Mori-Tanaka,  on  a  effectue  des  comparaisons 
avec  les  composites  de  2  et  3  phases.  Les  modeles  micromecaniques  ont  ete  evalues  avec  les 
donnees  de  comportement  mecaniques  disponsibles.  Les  comparaisons  ont  demontre  que 
le  modele  pour  les  vacuoles  a  calcule  correctement  les  comportements  pour  deux  types  de 
composites  charges  avec  0,2  jusqu’a  0,5  fraction  volume  de  particule.  Le  fait  qu’il  y  avait 
des  problemes  avec  les  predictions  de  module  initial  pour  le  composite  charge  avec  0,22 
fraction  volume  de  particule  indique  que  l’effet  de  la  non-linearite  de  la  matrice  devrait  etre 
inclus  dans  la  composition  actuelle. 
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EXECUTIVE  SUMMARY 


Propellants  are  presently  characterized  from  a  macroscopic  point  of  view.  This 
means  that  the  mechanisms  that  govern  material  behavior  are  lumped  together  and  mea¬ 
sured  as  a  unit  to  produce  a  material  property.  This  approach  does  not  provide  the  quan¬ 
titative  information  required  for  modifying  a  formulation.  To  meet  this  need,  an  analytical 
model  that  predicts  the  material  properties  from  knowledge  of  factors  such  as  particle  size 
distribution,  volume  fraction  of  particles,  adhesion  energy  and  polymer  properties  is  re¬ 
quired.  This  ability  to  predict  mechanical  properties  has  important  consequences  for  the 
determination  of  rocket  motor  service  life.  If  the  properties  of  the  motor  grain  can  be  pre¬ 
dicted  before  the  propellant  is  cast,  motor  service  life  can  be  determined.  If  the  calculated 
service  life  is  deemed  too  short,  the  model  can  be  used  to  guide  the  type  of  adjustments  that 
need  to  be  made  to  extend  the  service  life  of  the  motor.  This  capability  would  represent 
major  savings  in  development  and  life  cycle  management  costs  because  service  life  related 
problems  could  be  resolved  before  the  motor  is  fielded. 

In  recent  years,  researchers  in  the  propellant  industry  have  begun  to  use  composite 
materials  concepts  for  predicting  the  stress-strain  behavior  of  propellants.  These  concepts, 
based  on  a  microscopic  point  of  view,  take  into  account  the  size,  shape  and  quantity  of 
filler  introduced  into  polymeric  matrices.  Previously,  the  strengths  and  weaknesses  of  the 
Anderson-Farris  micromechanical  model  were  examined  using  experimentally  derived  data. 
This  model  was  based  on  an  energy  balance  concept  and  calculated  composite  modulus  using 
a  differential  scheme.  It  was  shown  that  the  model  could  predict  the  mechanical  behavior 
of  highly  loaded  composites  if  a  representative  value  of  adhesion  energy  was  available. 

In  this  report,  new  micromechanical  models  for  the  prediction  of  propellant  me¬ 
chanical  behavior  have  been  developed.  The  models  again  use  energy  balance  to  account 
for  nonlinear  behavior  due  to  particle  debonding.  However,  composite  modulus  is  predicted 
using  a  routine  based  on  an  improved  Mori-Tanaka  method.  This  method  permits  parti¬ 
cle  interaction  effects  to  be  taken  into  account  and  allows  the  stiffness  matrix  for  voids  or 
vacuoles  to  be  explicitly  stated  in  the  model.  Comparisons  with  literature  modulus  data 
showed  that  tensile  modulus  could  be  predicted  accurately  with  the  improved  Mori-Tanaka 
method.  Comparisons  with  literature  stress-strain  data  also  showed  that  the  model  derived 
for  vacuole  formation  predicted  the  mechanical  behavior  correctly  for  two  types  of  model 
propellants  which  contained  inclusion  volume  fractions  ranging  from  0.2  to  0.5. 
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NOMENCLATURE 

L  =  length  F  =  force 

square  brackets  [  ]  denote  dimensions  of  the  quantity 

A-ijki 

_ _ 

(Cijmn  ~  Gfjmn)-1  •  C°mnkl,  [ FL ~2] 

Bijkl 

“ 

(cg™»  -  CJwT1  •  [FL-*] 

c  ijki 

- 

average  elastic  constants  of  composite,  [FL  2] 

fo 

° ijkl 

- 

elastic  constants  of  comparison  material,  [FL~2] 

fir 

Uijkl 

- 

elastic  constants  of  phase-r  material,  [FL~2] 

E 

- 

volume  fraction  of  inclusions,  [— ] 

4 

- 

initial  volume  fraction  of  inclusions,  [— ] 

cr 

- 

volume  fraction  of  phase-r  inclusion,  [— ] 

cv 

- 

volume  fraction  of  voids  or  vacuoles,  [— ] 

E 

- 

average  composite  tensile  modulus,  [FL~2] 

Ec 

- 

average  composite  tensile  modulus,  [FL~2] 

Ei 

- 

inclusion  tensile  modulus,  [FL~2] 

Em 

- 

matrix  tensile  modulus,  [ FL~ 2] 

£ 

- 

composite  energy,  [FL] 

£int 

- 

interaction  energy,  [FL] 

£o 

- 

comparison  material  energy,  [FL] 

Go 

- 

adhesion  energy,  [FL] 

I ijkl 

- 

identity  matrix,  [— ] 

K 

- 

average  composite  bulk  modulus,  [FL~2] 

Pj 

- 

maximum  packing  fraction,  [— ] 

S 

- 

surface  area  of  inclusion,  [X2] 

or 

Tk 

ufj 

- 

Eshelby  matrix  of  phase-r  material,  [— ] 

- 

prescribed  surface  displacement,  [X] 

- 

constrained  displacement,  [X] 

V 

- 

volume  of  inclusion  [X3], 

Y 

- 

interaction  function  [— ], 

Ym 

- 

interaction  function  multiplier,  [— ] 

_ 

Kronecker  delta,  [— ] 

£cr 

- 

critical  strain,  [X/X] 

Gj 

- 

average  composite  strain,  [X/X] 

*ij 

- 

average  perturbed  strain,  [X/X] 
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e-  -  constrained  strain,  [L/L\ 

e*j  -  eigenstrain,  [L/L] 

e*j  -  corrected  eigenstrain  of  phase-r  material,  [L/L] 

e*Ju  -  uncorrected  eigenstrain  of  phase-r  material,  [L/L] 

r ijkl  -  correction  matrix  of  phase-r  material,  [— ] 

A ,/j,  -  Lame  constants,  [FL~2] 

/I  -  average  composite  shear  modulus,  [FL~2] 

V  -  average  composite  Poisson  ratio,  [— ] 

-  average  perturbed  stress,  [FL~2] 

a-j  -  prescribed  surface  stress,  [FL~2] 

afj  -  constrained  stress,  [FL~2] 

crfS  -  constrained  stress  of  phase-r  material,  [FL~2] 

-  inclusion  stress,  [FL~2] 

aljT  -  inclusion  stress  of  phase-r  material,  [Fir2] 


a*j  -  eigenstress,  [FL  2] 
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1.0  INTRODUCTION 

The  use  of  micromechanics  for  predicting  material  properties  is  not  new  in  the 
field  of  composite  materials  (Refs.  1-9).  In  recent  years,  researchers  in  the  solid  propellant 
industry  have  begun  using  composite  materials  concepts  for  predicting  the  stress-strain 
behavior  of  oxidizer /binder  type  propellants  to  assess  the  impact  that  a  formulation  may 
have  on  a  propellant’s  final  mechanical  properties  (Refs.  10  and  11). 

In  Refs.  12  and  13,  the  authors  hypothesized  that  nonlinear  mechanical  behavior  in 
filled  polymers  was  caused  by  the  debonding  of  inclusions  from  the  matrix.  They  modeled 
this  process  through  an  energy  balance  model  which  decreased  inclusion  concentration  and 
increased  void  concentration.  The  composite’s  mechanical  behavior  was  determined  by  the 
resulting  composite  modulus.  An  evaluation  of  their  technique  was  made  in  Refs.  14  and  15. 
It  was  concluded  that  their  model  could  predict  the  mechanical  behavior  of  highly  loaded 
composites  if  a  representative  adhesion  energy  was  available. 

Refs.  12  and  14  calculated  composite  modulus  using  the  Farber- Farris  (F-F)  routine 
(Ref.  16).  This  routine  is  an  incremental  scheme  based  on  a  model  that  was  developed  for 
the  prediction  of  viscosity  of  highly  concentrated  suspensions.  In  this  article,  the  F-F  routine 
has  been  replaced  by  a  modulus  prediction  routine  which  is  based  on  the  Mori- Tanaka 
(M-T)  method  (Ref.  17).  This  method  has  received  much  attention  recently  because  it 
permits  closed-form  solutions  for  multiphase  anisotropic  composites  (Refs.  18  and  19).  By 
including  recent  work  from  Ju  and  Chen  into  the  M-T  method,  particle-particle  interaction 
effects  can  be  taken  into  account  (Refs.  20  and  21).  This  results  in  an  improved  modulus 
prediction  routine  and  subsequently,  an  improved  micromechanical  model  for  the  prediction 
of  mechanical  behavior  of  particulate-filled  composites. 
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The  characteristics  of  the  improved  M-T  modulus  prediction  routine  will  be  exam¬ 
ined  and  compared  with  modulus  data  available  in  the  literature.  As  well,  comparisons 
will  be  made  with  F-F  predictions  to  illustrate  the  differences  between  the  two  routines. 
The  consequences  of  using  the  improved  M-T  routine  for  prediction  of  mechanical  behavior 
of  model  particulate  composites  will  be  shown  for  two  cases.  The  first  case  will  model 
debonded  particles  using  spherical  voids  to  simulate  completely  debonded  inclusions.  The 
second  case  will  model  debonded  particles  using  orthotropic  inclusion  properties  to  simu¬ 
late  spherical  inclusions  entrapped  by  spheroidal  voids.  Comparisons  will  be  made  between 
theoretical  predictions  of  mechanical  behavior  and  available  literature  data.  This  study  was 
undertaken  at  DREV  between  January  1994  and  March  1995  under  PSC  32C,  Rockets  and 
Missiles. 


2.0  COMPOSITE  MODULUS  PREDICTION 

Mori  and  Tanaka  showed  that  the  internal  stress  at  any  point  in  a  medium  con¬ 
taining  a  finite  concentration  of  inclusions  was  comprised  of  an  average  uniform  stress  and 
local  fluctuating  stresses  which  were  found  near  the  inclusions  (Ref.  17).  They  also  showed 
that  the  average  of  all  the  local  fluctuating  stresses  in  the  matrix  was  zero.  A  key  feature 
in  their  presentation  was  the  use  of  Eshelby’s  equivalent  inclusion  method  (Ref.  22). 

Given  the  importance  of  the  Eshelby  and  Mori- Tanaka  methods  to  the  current  work, 
it  would  be  appropriate  to  begin  with  an  overview  of  Eshelby’s  concepts  and  a  discussion 
of  the  salient  points  leading  to  the  Mori- Tanaka  formulation. 
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2.1  Eshelby’s  Equivalent  Inclusion  Method 


Eshelby  developed  a  theory  to  examine  how  a  far  field  uniform  applied  stress  or 
strain  is  disturbed  by  a  single  inclusion  in  an  infinite  elastic  matrix.  He  presented  it  in  the 
form  of  a  simple  series  of  cutting,  straining  and  rewelding  operations.  Conceptually,  the 
steps  were: 

1.  Remove  the  inclusion  from  the  matrix  and  allow  it  to  undergo  a  uniform  stress- free 
strain,  e*j.  The  stress-free  strain  is  also  known  as  the  transformation  strain  or  eigen- 
strain.  The  transformation  or  eigenstress  due  to  the  transformation  strain  would  be, 

=  A 4kSij  +  2  K- 

where  A  and  /i  are  the  matrix  Lame  constants. 

2.  Apply  a  surface  stress  —cr*jTij  to  the  inclusion  to  bring  it  back  to  its  original  shape. 
Embed  the  particle  back  into  the  matrix.  Consider  this  as  the  zero  strain  point  for 
the  inclusion  and  the  matrix. 

3.  Allow  the  surface  stress  —a*jUj  to  relax.  The  inclusion  will  see  a  change  in  stress  of 

+  AaijUj  while  the  matrix  will  see  an  equal  and  opposite  stress  of  ar*-n,j  — 

A CTijTij . 

The  fact  that  the  inclusion  had  different  properties  than  the  surrounding  matrix 
gave  rise  to  relaxed  or  constrained  strains,  e^-,  in  the  matrix  or  inclusion.  The  relaxed 
or  constrained  stresses  in  the  matrix  in  terms  of  the  matrix  properties  were  shown  to  be 
ofj  =  A €%hSij  +  2 iicij.  Thus,  the  stresses  in  the  inclusion,  a1- ,  were 


P501263.PDF  [Page:  16  of  122] 


UNCLASSIFIED 

4 


_C  _* 

aij  ~  aij 


Hekk  ~  €kk)dij  +  _  e*j) 


By  examining  the  elastic  field  in  an  ellipsoidal  inclusion,  Eshelby  arrived  at  an  important 
relationship  between  the  constrained  strains  and  the  eigenstrains, 


Cij  ~  Sijkl^kl 


where  Sijki  is  known  as  the  Eshelby  tensor.  The  tensor  is  a  function  of  the  matrix  Poisson 
ratio  and  the  shape  of  the  inclusion. 

Eshelby  had  also  shown  that  the  interaction  energy  between  the  constrained  elastic 
field  with  an  applied  elastic  field  uf  could  be  defined  as 


Sint  =  \Js(aiju?  ~  aijui)dSj 


where 


S  —  surface  area  of  inclusion, 

<7?.  =  constrained  inclusion  or  matrix  stress, 

erg  =  prescribed  surface  stresses, 

u'i  —  constrained  inclusion  or  matrix  strain, 

uf  =  prescribed  surface  displacements. 

Using  the  continuity  of  afj  and  ajj  across  the  inclusion  surface,  Gauss’  theorem,  the 
equivalence  and  eq.  1,  eq.  3  could  be  written  as 


Cm  = 

=  -\Sv^dv 


where 
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V 

e*- 


=  volume  of  inclusion, 

=  prescribed  surface  stresses, 
=  inclusion  eigenstresses, 

=  prescribed  surface  strains. 
=  inclusion  eigenstrains. 


The  interaction  energy  £,-nf  was  useful  because  it  could  be  used  to  find  the  elastic 
constants  of  a  composite.  Using  a  comparison  material  with  stiffness  properties  C°jkl~r 
along  with  eq.  5,  the  energy  of  the  composite,  £,  for  prescribed  stresses  was 


£  =  £0-£ 


int 


A  A 
aijGij 


^ /to  — 1  i 

2  ^ijkl  ai]  i 


rAdV 


[6] 


For  prescribed  strains,  the  energy  of  the  composite  was 


£ 


£o  +  £int 

\c?jk,44  -  I  Jv  r$4iV 


[7] 


The  meaning  of  the  energy  equations  may  be  interpreted  in  the  following  manner. 
For  hard  inclusions  (A,-  >  A0  and  m  >  fj,0),  the  inclusion  deforms  less  than  the  comparison 
material  for  a  prescribed  strain  ef-.  This  is  equivalent  to  Eshelby’s  case  where  an  inclusion 
spontaneously  contracts  (— e*j)  in  an  elastic  medium  due  to  phase  or  temperature  change. 
From  eq.  7,  the  negative  transformation  strain  signifies  that  the  composite  material  has 
an  increased  capacity  to  store  energy  in  relation  to  the  comparison  material.  Thus,  the 
elastic  constants  of  the  composite  are  greater  than  that  of  the  comparison  material.  For 
soft  inclusions  (A;  <  A0  and  /x,  <  /x0),  the  inclusion  deforms  more  than  the  comparison 
material.  Thus,  the  composite  has  a  decreased  capacity  to  store  energy.  This  leads  to  a 
composite  with  elastic  constants  which  are  less  than  those  of  the  comparison  material.  A 
similar  arguement  can  be  used  with  eq.  6. 
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Using  the  definitions  eki  =  Cijkl  1  cty  and  e°kl  —  C°jkl  eq.  6  may  be  used  to 
show  that  the  average  composite  strain  is 


€kl  —  4/  +  C€kl 


[8] 


where 

ekl  =  average  composite  strain, 

eh  =  comparison  material  strain, 

c  =  volume  fraction  of  the  inclusion, 

e*kl  —  eigenstrain  of  the  inclusion. 

Conversely,  it  can  be  shown  with  eq.  7  that  the  average  composite  stress  is 


Gij  =  Gij  -  caij 


where 

Wij  =  average  composite  stress, 

afj  =  comparison  material  stress, 

c  =  volume  fraction  of  the  inclusion, 

a*-  =  eigenstress  of  the  inclusion. 


[9] 


2.2  The  Mori- Tanaka  Method 


The  Mori-Tanaka  method  has  been  used  to  analyze  a  series  of  problems  ranging 
from  the  determination  of  composite  elastic  constants  (Ref.  19)  to  the  determination  of 
stresses  in  and  around  inclusions  (Ref.  23).  Following  Weng  (Ref.  18),  the  term  “phase” 
will  be  applied  to  a  collection  of  inclusions  whose  shape,  orientation  and  elastic  moduli  are 
identical.  The  superscript  r  will  be  used  to  identify  an  r-th  phase  of  inclusions. 

The  development  of  the  Mori-Tanaka  method  proceeds  as  follows  (Ref.  23).  For  a 
prescribed  stress  o-fi,  let  the  average  strain  produced  in  an  elastic  composite  and  comparison 
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material  be  defined  as 


„A 

aij 

—  C  ijkl^kl 

[10] 

„A 

aij 

II 

fro 

[11] 

where 

=  prescribed  surface  stresses, 

Cijkl 

=  average  elastic  constants  of  the  composite, 

tkl 

=  average  composite  strain, 

Cijkl 

=  elastic  constants  of  the  comparison  material, 

4i 

=  comparison  material  strain. 

Letting  the  average  perturbed  stress  due  to  the  presence  of  all  inclusions  be  defined 
as  &ij,  the  corresponding  perturbed  strains  are  ety.  The  overall  stress  in  the  comparison 
material  is 


atj  +  &ij  ~  C?jkl(€kl  +  *kl) 


[12] 


Looking  at  the  stresses  in  the  inclusion,  it  will  differ  from  the  stresses  in  the  compar¬ 
ison  material  by  an  amount  equivalent  to  Eshelby’s  constrained  stresses  due  to  the  difference 
in  material  properties.  The  inclusion  stresses  in  terms  of  the  inclusion  properties  are 

aij  —  afj  +  &ij  +  aij  =  +  Zkl  +  4i)  [13] 


Using  Eshelby’s  inclusion  method  eq.  1,  the  inclusion  stresses  in  terms  of  the  eigen- 
strains  are 


—  CTjj  +  &1] 


+  °ij  —  C°jkl(€kl  +  *kl  +  4l  —  €kl ) 


[14] 


Mori  and  Tanaka  showed  that  the  volume  integral  of  all  perturbed  stresses  in  a 
representative  volume  element  (RVE)  was  zero.  In  other  words,  the  perturbed  stresses  were 
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in  equilibrium  with  the  constrained  stresses.  Assuming  that  all  particles  in  the  RVE  were 
equally  stressed,  equilibrium  required 

°ij  +  cV.?/  =  0  [151 

Subtracting  eq.  12  from  eq.  14  and  using  eq.  2  for  phase  r, 

*ff  =  C!jkl(Sllmne*mnr-4n  [16] 

Defining  the  identity  tensor  as  I{jkl  =  +  &ilbjk)  and  using  eqs.  11,  12  and  16 

in  eq.  15,  the  perturbed  strains  are 

hi  =  -  £  „r  [IV] 

By  setting  the  third  term  in  eq.  13  equal  to  the  third  term  in  eq.  14  and  making  use 
of  the  definitions  for  comparison  strains  in  eq.  11,  constrained  strains  in  eq.  2  and  perturbed 
strains  in  eq.  17,  the  governing  equation  for  the  composite  is 

(c?jld-c;jtl)4i  =  (c&h  -  criHWklryr,r 

-  E  -  Atomic,']  +  cWt /  [is] 


2.3  Correction  for  Particle- Particle  Interaction 

Ju  and  Chen  (Ref.  21)  formulated  new  ensemble-averaged  constitutive  equations 
based  on  Eshelby’s  equivalent  inclusion  principle.  The  equations  took  into  account  particle- 
particle  interaction  effects.  Previously,  Weng  (Ref.  18)  had  shown  that  the  Mori-Tanaka 
(M-T)  method  could  be  reduced  to  give  the  Hashin-Shtrikman  (H-S)  lower  bound  as  a  spe¬ 
cial  case.  Experimental  data  has  shown  that  the  M-T  method  and  the  H-S  lower  bound 
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underestimated  slightly  composite  modulus  (Refs.  24  and  25).  Ju  and  Chen  showed  that 
by  accounting  for  particle  interaction,  their  model  could  reproduce  experimental  data  more 
closely.  As  a  special  case,  their  model  reduced  to  the  M-T  equations  when  particle  interac¬ 
tion  was  ignored.  This  meant  that  the  M-T  results  and  the  H-S  lower  bounds  corresponded 
to  a  micromechanical  model  with  non-interacting  particles. 


The  solution  in  Ref.  21  for  particle  interaction  comes  in  a  form  convenient  for 
inclusion  in  a  M-T  formulation.  Since  the  theory  is  based  on  Eshelby’s  equivalent  inclusion 
principle,  the  correction  matrix,  F",  could  be  cast  in  terms  of  the  un corrected  eigenstrain 
solution.  From  this  point  on,  braces  will  be  used  to  denote  vectors  and  brackets  will  be 
used  to  denote  square  matrices.  The  corrected  eigenstrain  vector  of  phase-r  was  shown  to 
be 

ter}  =  [rrKd  [is] 

where 

Rr)  = 

in 

Rr>  = 

[i] 

cr  = 

Y 

[FT] 

Ci  = 

C2  = 

a  = 

P 

VG  = 

K0,Kr  = 

f^oi  f^r  — 


interacting  or  corrected  eigenstrain  of  phase-r, 

m+f fY[w% 

non-interacting  or  uncorrected  eigenstrain  of  phase-r, 

identity  matrix, 

volume  fraction  of  phase-r, 

interaction  factor, 

Cl  Gijtikl  +  CiihkSjl  +  SuSjk), 

12(13z/c  -  14u2)  -  ^(1  -  2va)(l  +  v0), 

6(25  -  34u0  +  22^)  -  3^,(1  -  2v0)(l  +  u0), 

2(5^o  —  1)  +  10(1  -  Vo)  Ik^-Ko  ~  ’ 

2(4  -5^)+  15(1  -v0){-^), 

Poisson’s  ratio  of  comparison  material, 

bulk  modulus  of  comparison  and  phase-r  inclusion, 

shear  modulus  of  comparison  and  phase-r  inclusion. 


In  eq.  19  the  inter-particle  effects  are  quantified  by  the  second  term  in  pT].  As  ex¬ 
plained  in  Ref.  21,  this  matrix  was  derived  from  the  analysis  of  probabilistic  pairwise  particle 
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interaction  of  two  identical  and  randomly  located  elastic  spheres  which  are  embedded  in  a 
comparison  material. 

The  factor  Y  (Ref.  21)  is  related  to  the  microstructure  of  the  composite  and  is 
characterized  by  a  radial  distribution  function,  g(r).  For  a  statistically  uniform  radial 
distribution  function  where  g(r )  =  1,  Y  =  1/24.  It  will  be  shown  later  that  |Tr]  causes 
certain  problems  in  the  recovery  of  [ CT ]  when  cr  =  1.  This  can  be  remedied  by  selecting  a 
modified  form  of  Y.  This  point  will  be  discussed  further  in  Sec.  6.0. 


TWO-PHASE  COMPOSITE 


The  average  composite  modulus  for  an  elastic  medium  containing  a  finite  concen¬ 
tration  of  inclusions  may  be  calculated  using  the  energy  equation  that  Eshelby  developed 
but  extended  for  multiple  inclusions.  If  the  domain  of  integration  in  eq.  6  is  extended  for 
finite  concentrations  of  inclusions  of  phase  r,  then  the  following  equation  can  be  written 


A  A  _ 
aklaij  ~ 


r<o 

2  '“'ijkl 


-e*rdVT 


This  leads  naturally  to  a  modified  form  of  eq.  8  which  denotes  that  all  particles  of  phase  r 
are  taken  into  account, 

«  =  K}  +  X>r{£H  [2i] 

where  {e*r}  is  defined  in  Sec.  2.3.  The  form  of  eq.  21  is  slightly  different  than  that  found 
in  Ref.  18.  Here  the  interacting  eigenstrain  is  used.  However,  if  |Tr]  is  set  equal  to  [I],  the 
non-interacting  eigenstrain  is  easily  recovered.  This  would  give  an  average  composite  strain 
identical  to  that  found  in  Ref.  18. 


The  average  elastic  constants  for  a  two-phase  composite  can  be  determined  in  a 
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two-step  process  using  the  equations  discussed  in  Secs.  2.2  and  2.3.  First,  eqs.  18  and  19 
are  used  in  eq.  21  to  obtain 

{e*r}  =  _(cr[J-5r]  +  [Sr]  +  {Cr-C°}-1  •  [ C° ]  -  cr[r'])"1{<f}  [22] 

Note  that  c^r  in  eq.  18  is  identical  to  {e*r}  in  eq.  19.  Second,  solving  eq.  21  for  {e0}  and 
substituting  it  into  eq.  11  gives 

{”*}  =  [C°]  •  ({f}  -  c'K’-})  [23) 

Then  using  eqs.  19  and  22  into  eq.  23  and  making  a  sign  change,  the  average  elastic  properties 
can  be  expressed  as 

[C]  =  [C°]  •  ([/]  +  cr[rr](cr[J-Sr-rr]  +  [5r]  +  [C-C0]-1  •  [C0])"1)  [24] 

Comparison  of  eq.  24  with  that  of  eq.  53  in  Ref.  21  shows  that  the  two  are  not 
identical  for  the  case  where  particle  interaction  is  taken  into  account.  The  difference  occurs 
because  Ref.  20  derived  the  volume- averaged  strain  tensor  by  integrating  Green’s  function 
over  a  representative  volume  element  to  arrive  at  a  definition  that  was  a  function  of  the 
uniform  far-held  strains,  the  eigenstrains  and  a  depolarization  tensor. 

4.0  ELASTIC  PROPERTIES  OF  A  THREE-PHASE  COMPOSITE 

In  order  to  predict  the  mechanical  behavior  of  a  particulate  composite,  the  model 
in  Refs.  12  and  14  required  new  elastic  moduli  be  calculated  each  time  a  group  of  inclusions 
are  debonded.  The  elastic  modulus  of  a  composite  with  bonded  and  debonded  inclusions 
was  assumed  equivalent  to  the  elastic  modulus  of  a  composite  which  contains  inclusions 
and  voids.  Thus,  the  case  of  a  three-phase  composite  is  of  special  interest  here  because  the 
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deboriding  process  in  the  above  assumption  changes  the  initial  two-phase  composite  where 
all  inclusions  are  bonded  to  the  matrix  into  a  three-phase  composite  where  some  particles 
become  debonded  thereby  creating  voids  or  vacuoles.  The  distinction  between  a  void  and 
vacuole  is  that  a  void  is  a  spherical  air  bubble  while  a  vacuole  is  a  spheroidal  air  pocket 
which  surrounds  a  debonded  inclusion. 

The  average  elastic  properties  for  a  3- phase  composite  can  be  derived  using  the 
procedure  outlined  in  Sec.  3.0.  The  only  difference  is  that  the  summation  of  phases  will 
now  be  used.  Therefore,  eq.  18  is  modified  by  explicitly  specifying  the  number  of  phases  in 
eq.  17.  The  two  resulting  equations  from  consideration  of  each  phase  are 

-{e°}  =  (c*[I— 5‘]  +  [S'4']  +  [A]){e*1'}  +  cv[I-Sv]{e*v}  [25] 

~{e°}  =  (cu[7-5v]  +  [5,;]-l-[JB]){c*u}  +  ct'[J-5i]{e**'}  [26] 

where 

[A]  =  [Ci  —  C°]~1  •  [C°], 

[5]  =  [Cv-C0}~1  -[C% 

i  =  parameters  relating  to  inclusions, 

v  =  parameters  relating  to  voids  or  vacuoles. 

Subtracting  eq.  26  from  eq.  25  gives  the  relationship  between  {e**}  and  {e*“}, 

([S’*']  +  [A]){e**'}  -  ([SI  +  [B]){e"’}  -  {0}  [27] 

Following  the  procedure  outlined  in  Sec.  3.0,  the  non-interacting  eigenstrains  for 
the  inclusion  and  void  or  vacuole  are 

{<?'}  =  -(cii-s4'-r4']  +  [s‘']  +  [A] 

+cv[I-Sv-Tv]  •  [£v  +  B]~l  •  [S{  +  A])-1{e} 


[28] 
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(C>  =  -(cv[I-Sv-Tv]  +  [Sv]  +  [B] 

+ci[I-Si-Ti ]  •  [5*'  +  A]"1  •  [Sv  +  -B])-1{e}  [29] 

As  before,  starting  from  eq.  11,  the  substitution  of  eqs.  21,  19,  28  and  29  gives  the 
average  elastic  properties  of  a  3-phase  composite  as 

[c]  =  [c0]  •  ([/]  +  ct[r](ci[/-s”'-ri]  +  [si  +  [A] 

+cv[I-Sv-Tv]  •  [S*  +  B l"1  •  [S'*  +  A])"1 
+  c'J[ri(ci7-S',;-ri  +  [5i  +  [JB]  [30] 

+ci[/-5i-r]  •  [S{  +  A]-1  •  [S'"  +  B})-1) 

As  mentioned  in  Sec.  2.1,  the  Eshelby  tensor  [5]  is  dependent  on  the  matrix  Poisson 
ratio  u0  and  the  inclusion  shape.  Reference  23  gives  [5]  for  spheroidal  inclusions.  For 
spherical  inclusions  [5]  is  defined  by  Ref.  21  as 

Sijki  =  15(1  _  v  ~  1  )&ijf>ki  +  (4  -  5v0)(6ikSji  +  SuSjk ))  [31] 

For  the  case  where  the  debonded  inclusions  are  modeled  by  equivalent  sized  spherical 
voids,  the  property  matrix  [ Cv ]  can  be  set  to  zero.  Alternatively,  isotropic  void  properties 
can  be  specified  to  simulate  partially  debonded  inclusions.  This  technique  was  used  in 
Ref.  12.  In  either  case,  the  average  composite  properties  remains  isotropic  so  the  usual 
engineering  elastic  constants  like  bulk,  shear  and  tensile  modulus  and  Poisson’s  ratio  may 
be  found  from  eq.  30  using 
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K  =  (l/3)(C'iin  +  2(7x122) 

F  =  (l/2)(C'nii  +  (7ii22) 

E  =  E„„  -  [32] 

L'  2222  +  (^2233 
_  _  (72211 

^2222  +  C 2233 

In  the  case  of  modeling  vacuoles,  Mochida  (Ref.  26)  suggested  that  the  boundary 
conditions  around  a  debonded  inclusion  may  be  specified  as 

0  =  0%  =  Ciiki(e°kl  +  hi  +  (%i)  [33] 

0  =  €jj  +  €jj  +  €jj  [34] 

Assuming  the  load  is  applied  in  the  n-direction  and  the  equator  of  the  particle  is  represented 
by  the  jj-direction,  eqs.  33  and  34  state  that  no  stresses  are  experienced  by  the  inclusion 
in  the  loading  direction  and  no  lateral  contraction  is  allowed  around  the  equator. 

Attempting  to  apply  Mochida’s  boundary  conditions  to  the  present  formulation 
results  in  a  trivial  solution  where  {e*}  =  {0}  because  application  of  eqs.  33  and  34  to 
eqs.  13  and  14  leads  one  to  conclude  that 

{0}  =  {e°  +  €  +  h} 

{0}  =  {€°  +  z  +  ec  -  e*} 

An  alternative  to  Mochida’s  method  is  to  model  the  debonded  phase  as  a  spherical 
inclusion  with  orthotropic  properties.  A  low  or  zero  modulus  value  in  the  loading  direction 
can  be  used  to  represent  the  debonded  condition  and  a  high  or  inclusion  modulus  value 
in  the  equator  direction  can  be  used  to  enforce  the  lateral  constraint  condition.  Since  the 
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M-T  formulation  can  be  applied  equally  well  to  inclusions  with  orthotropic  properties  as 
to  inclusions  with  isotropic  properties,  this  approach  can  be  implemented  by  modifying  the 
definition  of  the  debonded  particle’s  material  matrix.  The  property  matrix  for  the  normal 


components  of  an  orthotropic  material  is  (Ref.  27) 


where 


[Cv]  =  m 


E\ i(l  -  ^23^32) 
£22(^12  +  ^13^32) 
£33(^13+^12^23) 


£11(^21  +^23  ^3l) 
£22(1  —  ^13^31) 
£33(^23  +  ^13^21) 


£ll(^31  +  ^21^32) 
£22(^32  +  ^12^31) 
£33(1  —  ^12^21) 


m  =  (1  —  l/12t'21  ~  ^13^31  —  ^23^32  —  2/12?/23^/31  ~  ^13^21  ^32 )_1 1 

Eu  =  tensile  modulus  of  vacuole  in  the  u-direction, 

Uij  =  Poisson’s  ratio  of  vacuole  in  the  indirection. 


[35] 


5.0  CRITICAL  STRAIN  FOR  ORTHOTROPIC  COMPOSITES 


The  prediction  of  mechanical  behavior  is  performed  using  an  energy  balance  model 
which  is  derived  in  terms  of  critical  strain.  Critical  strain  is  defined  as  the  point  where  the 
internal  strain  energy  in  the  composite  and  the  energy  released  due  to  particle  debonding 
equals  the  work  put  into  the  composite.  This  statement  can  be  expressed  as  (Refs.  12  and 
14) 

2  GC6A!V0  =  crijSeij  —  6aij€ij  [36] 

where  Gc  is  the  adhesion  energy  between  particle  and  matrix,  8  A  is  the  variation  in  surface 
area,  crtj  is  the  composite  stress,  Cjj  is  the  composite  strain  and  V0  is  the  specimen  volume. 

Examining  the  case  of  a  uniaxial  bar  under  tension  and  ambient  pressure,  eq.  36 
can  be  simplified  to  (Refs.  12  and  14) 


2 GC8A/V0  =  a^Scer  - 


[37] 
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where  is  the  applied  stress  in  the  11-direction  and  ecr  the  critical  strain.  The  boundary 
conditions  <722  =  033  —  0  and  €22  =  e33  were  used. 

By  using  the  same  boundary  conditions  on  the  constitutive  equation  for  an  or¬ 
thotropic  material  (eq.  10),  the  transverse  strain  is 

C2211 


e22  =  -7=7 


C 2222  +  C 2233 

Using  eq.  38  to  solve  for  aii  in  eq.  10,  the  simplified  constitutive  equation  is 


[38] 


&11  —  C’lin  — 


2C1122C2 


-= -  €, 


C 2222  +  C 2233, 

Taking  the  variation  of  cr^  and  substituting  it  and  eq.  39  into  eq.  37  gives 

f  Ci  122^221 1(^2222  +  8C  2233)  1 


[39] 


2  GC6A/V0  = 


-8C1111  -  2 


7 


[40] 


9  J  (C2222  +  C,2233)(Cr2211^Cn22  +  1122$ C 221l)  1 

1  /7=f  .  ~  \0  |  I 


(C 2222  +  C2233)2 

S3)(Cr2211^Cn22  + 

(^2222  +  C2233)2 

This  equation  is  valid  for  orthotropic  composites  under  uniaxial  tension  and  ambient  pres¬ 
sure.  It  assumes  the  RVE  is  larger  than  the  largest  particle  so  that  average  stress,  strain 
and  modulus  can  be  used.  If  [C]  is  isotropic  then  eq.  40  reduces  to 


2GcbAjV0  =  -6Ee, 


[41] 


where  6E  is  the  variation  in  tensile  modulus. 

To  solve  eq.  40  or  eq.  41,  the  variational  quantities  are  changed  to  incremental 
quantities  based  on  inclusion  concentration.  For  example,  6  A  would  change  to  A  A/  Ac*.  The 
algorithm  to  generate  the  composite’s  stress-strain  behavior  has  been  outlined  in  Refs.  14 


and  15. 
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6.0  MODULUS  PREDICTION  PERFORMANCE 

The  performance  of  the  energy  balance  model  in  Ref.  14  is  directly  related  to  the 
performance  of  the  composite  modulus  prediction  routine  employed.  Therefore,  it  is  of  in¬ 
terest  to  examine  the  prediction  characteristics  of  the  composite  modulus  models  developed 
in  Secs.  3.0  and  4.0.  To  accomplish  this  task,  the  models  will  be  compared  with  experi¬ 
mental  data  obtained  from  Refs.  25,  29  and  30.  As  well,  they  will  be  compared  with  the 
Farber-Farris  predictions  (Ref.  16).  The  Farber-Farris  routine  was  originally  used  in  the 
energy  balance  model  in  Ref.  14.  The  Mathematica  (Ref.  28)  input  file  used  for  genera¬ 
tion  of  the  Mori- Tanaka  data  may  be  found  in  Appendix  A.  A  summary  of  the  constituent 
properties  for  the  composites  used  in  this  document  are  shown  in  Table  I. 

6.1  Comparison  with  Glass/Epoxy  Composite 


A  comparison  of  the  2-phase  composite  predictions  (eq.  24)  with  experimental  data 
obtained  from  Smith  (Ref.  29)  and  theoretical  predictions  using  the  Farber-Farris  routine 
is  shown  in  Fig.  1.  Smith’s  composite  system  consisted  of  glass  spheres  embedded  in  epoxy 
( Ei/Em  =  25).  Three  variations  of  the  2-phase  model  are  shown. 

The  first  variation  is  the  Mori- Tanaka  model  where  no  particle  interaction  is  ac¬ 
counted  for  (denoted  M-T  in  Fig.  1,  [rr]  =  [/]  in  eq.  24).  This  variation  will  be  called 
simply  the  M-T  model.  It  can  be  seen  that  the  M-T  model  underpredicted  Smith’s  data  at 
volume  fractions  c*  above  0.3.  This  was  not  surprising  since  Ref.  18  showed  that  the  M-T 
solution  was  equivalent  to  the  Hashin-Shtrikman  lower  bounds. 


The  second  variation  (called  M-T  interaction)  is  the  M-T  model  corrected  for  par- 
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tide  interaction  ([Tr]  with.  Y  —  1/24).  This  variation  fared  better  up  to  c*  =  0.4.  However, 
between  0.4  <  c*  <  0.63,  the  normalized  composite  modulus,  Ecj Em ,  rapidly  rose.  At 
c*  =  0.63,  it  dropped  suddenly.  Between  0.63  <  c*  Cl,  Ec/Em  rose  back  up  to  —2.  It 
was  obvious  that  this  undesirable  behavior  was  being  caused  by  the  introduction  of  the 
correction  matrix  [rr].  An  examination  of  eq.  24  showed  that  unless  [rr]  —»•[/]  as  c*  —*■  1, 
the  inclusion  properties  could  never  be  recovered.  One  way  of  forcing  this  behavior  was  to 
introduce  a  factor  1  —  cl  into  the  definition  of  Y  in  eq.  19. 

Thus,  the  third  variation  (called  M-T  modified)  is  the  M-T  model  corrected  using  a 
modified  partide  interaction  factor  ([rr]  with  Y  =  (1  —  c*)/ 24).  This  modd  behaved  much 
better  and  gave  good  predictions  up  to  c*  =  0.5.  In  comparison  with  the  Farber-Farris 
predictions  in  which  a  maximum  packing  fraction  Pf  =  1  was  used  (denoted  F-F  in  Fig.  1), 
the  M-T  modified  model  gave  similar  results. 

An  attempt  was  made  to  see  if  Y  =  (l-c*)/24  had  any  significance  when  compared 
with  the  traditional  radial  distribution  function  predicted  by  the  Percus-Yevick  solution 
(Refs.  31  and  32).  A  comparison  showed  there  was  no  similarity  between  the  two.  Thus, 
the  function  Y  used  here  only  provides  a  means  of  recovering  inclusion  properties  and  cannot 
be  seen  as  an  indicator  of  microstructural  features  as  suggested  by  the  authors  in  Ref.  21. 

6.2  Comparison  with  Tungsten-Carbide/Cobalt  Composite 


A  similar  comparison  is  made  using  experimental  data  obtained  from  Ravichandran 
(Ref.  25).  Here  the  composite  system  is  a  tungsten-carbide/cobalt  cermet  (Ei/Em= 3.4). 
The  results  are  shown  in  Fig.  2. 
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In  contrast  to  Fig.  1,  Fig.  2  shows  that  the  predicted  j Ec/ Em  rises  smoothly  for 
the  entire  range  of  c*.  This  is  a  consequence  of  the  low  inclusion  to  matrix  modular  ratio. 
This  material  system  demonstrates  the  sensitivity  of  the  M-T  interaction  solution  to  the 
constituent  properties. 

The  M-T  modified  and  Farber-Farris  solutions  give  almost  identical  results  while 
the  M-T  solution  gives  slightly  lower  modulus  values.  These  three  models  tended  to  un¬ 
derpredict  the  experimental  data.  The  M-T  interaction  solution  shows  that  it  overpredicts 
modulus  in  the  neighborhood  of  c'  =  0.75  and  above. 

6.3  Comparison  with  Glass/Polyurethane  Composite 


As  a  final  comparison,  the  experimental  data  from  Yilmazer  is  used  (Ref.  30).  The 
glass  bead/polyurethane  composite  had  an  JS,-/ Em  >  16000.  This  material  system  exhibited 
behavior  that  was  different  than  the  previous  two.  At  cl  =  0.5,  the  normalized  modulus 
measured  by  Yilmazer  was  Ec/ Em  =  13.  When  other  references  (Ref.  25  in  Ref.  25,  Refs. 
29,  33,  34  and  35)  were  examined,  generally  Ec/ Em  <  5  at  c*  =  0.5. 

The  Farber-Farris  solution  underpredicted  the  experimental  results  when  the  max¬ 
imum  packing  fraction  Pf  —  1  was  used  (see  Fig.  3).  To  bring  the  F-F  solution  up  to 
the  measured  data,  Pj  was  set  equal  to  0.6.  The  fact  that  a  packing  fraction  is  required 
indicates  that  the  additional  reinforcement  seen  in  this  material  system  is  not  due  to  the 
Ei/Em  alone  but  includes  some  other  reinforcing  mechanism  like  particle  interaction. 

As  expected  the  M-T  solution  underpredicted  the  experimental  data.  Again,  to 
reproduce  the  experimental  data,  a  modification  to  the  definition  of  Y  in  eq.  19  was  required. 
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This  time,  an  interaction  factor  multiplier  Ym  =  1.26  was  added  giving  Y  —  Ym(  1  —  c*)/ 24. 
The  value  for  Ym  was  found  through  trial  and  error.  Like  the  maximum  packing  fraction, 
Ym  provides  an  indication  of  additional  reinforcement  through  particle  interaction.  As  can 
be  seen  from  Fig.  3,  the  M-T  modified  solution  reproduces  the  measured  moduli  quite  well. 

6.4  Comparison  for  Three-Phase  Composites 

Experimental  data  for  the  modulus  of  three-phase  composites  is  scarce  (Ref.  18). 
A  study  by  Ishai  and  Cohen  on  a  composite  system  comprised  of  sand  and  voids  in  a  matrix 
of  epoxy  is  frequently  used  for  evaluation  of  three-phase  models  (Ref.  35).  The  modular 
ratio  for  this  system  is  Ei/Em  =  36.  Huang  (Ref.  36)  showed  that  the  volume  fraction  of 
sand  c*  varied  with  volume  fraction  of  voids  cv  according  to  c*  =  0.173(1  —  cv)  for  the  Ishai 
and  Cohen  data.  With  composites  containing  void  volume  fractions  of  0.1  <  cv  <  0.6,  the 
inclusion  volume  fraction  was  within  a  narrow  range  of  0.16  >  c*  >  0.07. 

Figure  4  shows  a  comparison  of  the  Farber-Farris  (F-F)  model  and  the  3-phase  M-T 
model  (eq.  30).  Two  variations  of  the  M-T  model  have  been  used.  The  first  is  the  M-T 
model  with  no  particle  interaction  (denoted  M-T  in  Fig.  4).  The  second  is  the  modified 
M-T  model  with  particle  interaction  between  inclusions  and  between  voids  (denoted  M-T, 
P,  T").  An  interaction  factor  of  Y  =  (1  -  c‘)/24  was  used  for  the  inclusions  and  voids. 
As  Fig.  4  shows,  all  models  predicted  the  experimental  composite  moduli  well.  This  was 
expected  since  the  inclusion  loading  was  fairly  low.  Small  differences  between  the  three 
models  can  be  seen  though.  The  F-F  model  tended  to  predict  higher  moduli  at  the  lower 
void  concentrations  and  vice  versa  at  the  higher  void  concentrations.  The  M-T  model 
predicted  higher  moduli  than  the  F-F  model  and  the  modified  M-T  model  at  the  higher 


P501263.PDF  [Page:  33  of  122] 

UNCLASSIFIED 

21 

void  concentrations.  The  modified  M-T  model  predicted  the  experimental  data  over  the 
entire  range  except  at  cv  =  0.4  where  it  overpredicted  the  composite  modulus. 

A  further  understanding  of  the  3-phase  M-T  model  is  gained  by  predicting  composite 
modulus  using  the  Yilmazer  glass  bead/polyurethane  (Ref.  30)  and  Smith  glass  bead/epoxy 
(Ref.  29)  property  data.  Since  the  energy  balance  model  (Sec.  5.0)  assumes  debonded 
inclusions  will  become  voids  or  vacuoles,  the  void  or  vacuole  volume  fraction  varies  according 
to  cv  =  c‘D  —  c*.  An  initial  inclusion  volume  fraction  c*  of  0.5  has  been  selected  here.  Voids 
were  modeled  by  setting  [( Cv )  =  [0].  Vacuoles  were  modeled  by  setting  =  0  and  i/”2  = 
z/J'g  =  0  in  eq.  35  assuming  loading  was  in  the  11-direction.  A  value  of  Y  =  1.26(1  —  c‘)/24 
was  used  for  the  generation  of  the  Yilmazer  M-T  void  and  M-T  vacuole  predictions  while 
a  value  of  Y  =  (1  -  c*)/ 24  was  used  for  the  Smith  predictions.  The  interaction  multipliers 
were  chosen  based  on  the  results  obtained  in  Sec.  6.3. 

Figure  5  shows  the  performance  of  the  F-F  and  M-T  models  for  the  glass/void/ 
polyurethane  (PU)  data.  The  M-T  solution  (no  interaction)  is  shown  for  reference.  For 
both  the  M-T  void  and  M-T  vacuole  results,  a  larger  drop  in  modulus  is  predicted  when 
compared  to  the  F-F  results  (Pj  =  0.6).  However,  at  c1  <  0.17,  all  models  gave  similar 
moduli.  The  M-T  vacuole  model  gave  slightly  higher  moduli  than  the  M-T  void  model  as 
c‘  — >0  (see  Evac/ Evoid  curve).  This  was  expected  since  the  vacuole  stiffness  matrix  [Cv\  was 
partially  populated  with  non-zero  values.  This  behavior  was  also  observed  by  Mochida  in 
Ref.  26. 

The  model  performance  for  the  glass/void/epoxy  data  (Fig.  6)  showed  that  as  c‘  — ►  0, 
the  M-T  models  predicted  higher  modulus  between  0.45  <  c‘  <  0.5  then  lower  moduli  down 
to  cl  —  0.12  in  comparison  to  the  Farber- Farris  solution  ( Pf  =  1).  This  was  not  surprising 
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given  the  2-phase  predictions  in  Sec.  6.1.  The  M-T  model  moduli  decreased  more  rapidly 
than  the  F-F  modulus  but  the  difference  in  rates  of  decrease  was  not  as  great  as  that  seen 
for  the  glass/PU  data.  Again,  the  M-T  vacuole  model  gave  a  slightly  higher  moduli  than 
the  M-T  void  model  as  cl  — >  0. 

A  comparison  of  predicted  composite  Poisson  ratio  vc  for  the  glass/void/PU  com¬ 
posite  is  shown  in  Fig.  7.  The  F-F  predictions  gave  a  slightly  concave-down  curve  as  c*  — »-0. 
The  M-T  solutions  on  the  other  hand  had  a  curious  concave-up  shape. 

The  dependency  of  the  Poisson  ratio  curvature  on  the  matrix  Poisson  ratio  vm  is 
shown  in  Fig.  8.  The  solid  lines  represent  the  Poisson  behavior  if  the  inclusions  were 
assumed  to  be  rigid  {[&)  —  [oo]).  Reference  20  showed  that  vc  was  independent  of  the 
inclusion  Poisson  ratio  for  this  case.  For  a  matrix  Poisson  ratio  vm  =  0.2,  it  can  be  seen 
that  the  composite  Poisson  ratio  vc  equalled  0.2  for  all  cJ.  Budiansky  (Ref.  37)  observed 
the  same  behavior  in  his  analysis  of  elastic  moduli  of  heterogeneous  materials.  He  found 
that  at  c4  =  0.2,  his  shear  and  bulk  modulus  equations  became  decoupled.  Calculations  for 
I'm  >  0.2  show  that  with  rigid  inclusions,  uc  has  a  concave-up  shape.  For  um  <  0.2,  vc  has  a 
concave-down  shape.  It  appears,  then,  for  analyses  which  use  a  discrete  averaging  process 
such  as  that  given  by  eqs.  24  and  30  or  employed  in  Refs.  18,  20  or  37,  the  predicted  vc 
will  show  this  kind  of  behavior. 

The  dashed  lines  in  Fig.  8  show  the  vc  behavior  using  glass/PU  property  data.  vm 
was  varied  between  0.1  to  0.3  while  i \  =  0.16.  Particle  interaction  was  accounted  for  and 
Y  —  1.26(1  —  c')/ 24.  The  deviation  from  the  3-phase  rigid  inclusion  prediction  at  vm  =  0.2 
is  caused  by  v{  /  0.2.  For  values  of  vm  ^  0.2,  the  effects  of  |Tr]  and  Ei/Em  come  into  play. 
Since  vm  =  0.499  for  PU,  it  is  evident  why  a  concave-up  curve  is  predicted  in  Fig.  7. 
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Going  back  to  Fig.  7,  the  M-T  solution  gives  a  higher  prediction  of  vc  compared  with 
the  M-T  void  and  M-T  vacuole  solutions.  In  the  case  of  the  M-T  void  results,  vc  eventually 
converges  to  the  M-T  solution  since  c*  — *•  0.  The  M-T  vacuole  results  demonstrate  the 
restraining  power  of  the  —  v\z  —  0  assumption.  Numerical  trials  for  glass/ void/epoxy 
showed  the  same  trends. 

7.0  PREDICTION  OF  MECHANICAL  BEHAVIOR 

The  general  characteristics  of  the  critical  strain  equation  (eq.  40)  will  first  be 
examined  using  the  M-T  void  model.  Stress-strain  behavior  for  different  volume  fractions 
of  inclusions  will  be  shown  to  illustrate  how  the  model  functions.  Void  volume  fraction  is 
calculated  from  cv  —  cl0  —  c*  where  cx0  is  the  initial  inclusion  volume  fraction. 

To  examine  the  consequences  of  using  the  critical  strain  equation  based  on  the 
M-T  void  or  M-T  vacuole  models,  mechanical  behavior  predictions  will  be  compared  with 
literature  data  for  two  types  of  model  particulate  composites.  The  first  composite  is  a 
glass  bead/polyurethane  system  which  was  studied  by  Yilmazer  and  Farris  (Ref.  30).  The 
second  composite  is  a  glass  bead/polyethylene  system  which  was  studied  previously  by  the 
author  (Ref.  14).  Predictions  using  the  Farber-Farris  model  will  be  made  to  highlight  the 
differences  between  it  and  the  M-T  models.  The  FORTRAN  program  used  for  generating 
the  M-T  void  and  vacuole  data  may  be  found  in  Appendix  B. 
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7.1  Model  Characteristics 

The  results  shown  in  Figs.  9  and  10  were  generated  using  the  constituent  properties 
from  Ref.  29  (see  Table  I).  Y  =  (1  -  c*)f 24  was  used  (Sec.  6.1)  along  with  an  adhesion 
energy  of  5  J/m^  (Ref.  13).  Particle  diameter  size  was  set  at  50  /im  and  the  log  standard 
deviation  equalled  0.2. 

Figure  9  shows  how  stress-strain  behavior  changes  with  initial  inclusion  volume  frac¬ 
tion.  As  initial  volume  fraction  increased,  composite  modulus  also  increased.  This  trend 
has  been  well  documented  in  the  literature  (Refs.  38,  39  and  40).  Examining  the  volume 
fraction  curves,  as  long  as  initial  volume  fraction  remained  constant,  linear  elastic  behavior 
was  predicted.  Maximum  stress  was  determined  by  the  composite  modulus  and  the  initial 
critical  strain  values.  Once  inclusions  started  to  debond,  volume  fraction  decreased  and 
nonlinear  behavior  was  observed.  Stress  increased  or  decreased  depending  on  the  initial 
volume  of  inclusions.  Below  inclusion  volume  fractions  of  0.05,  linear  behavior  was  again 
predicted.  The  modulus  at  this  point  was  determined  by  the  volume  of  voids  in  the  compos¬ 
ite.  The  current  model  does  not  contain  a  failure  criterion  so  only  a  portion  of  the  predicted 
stress-strain  curve  would  be  seen  in  experiments. 

The  relationship  between  volume  fraction  and  probability  of  inclusion  survival  with 
strain  is  shown  in  Fig.  10.  Probability  of  survival  represents  the  number  of  well-bonded 
inclusions  remaining  in  the  composite.  It  can  also  be  thought  of  as  a  cumulative  size 
distribution  curve  where  the  lefthand  side  of  the  figure  represents  large  particle  diameters 
and  the  righthand  side  representing  small  diameters.  For  the  log  standard  deviation  studied, 
10%  of  the  total  number  of  particles  made  up  50%  of  the  inclusion  volume.  When  the 
inclusion  volume  fell  to  <  0.05,  approximately  50%  of  the  total  number  of  particles  remained. 
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This  shows  that  the  larger  particles  controlled  the  type  of  nonlinear  behavior  predicted  and 
that  the  smaller  particles  were  almost  inconsequential. 

7.2  Glass  Bead/Polyurethane  Composite 


In  Ref.  30,  the  authors  fabricated  glass  bead/polyurethane  composites  with  in¬ 
clusion  volume  fractions  ranging  from  c\  —  0.0  to  cl0  —  0.5  in  steps  of  0.1.  One  set  of 
composites  contained  untreated  glass  beads  while  a  second  set  used  glass  beads  treated 
with  a  silane  coupling  agent.  Six  out  of  the  twelve  composite  combinations  studied  in  Ref. 
30  have  been  chosen  to  demonstrate  the  performance  of  the  models.  In  Ref.  13,  the  authors 
used  void  properties  and  adjusted  the  matrix  modulus  and  the  adhesion  energy  to  obtain 
agreement  between  the  experimental  and  predicted  values.  In  this  section,  only  adhesion 
energy  was  considered  an  adjustable  parameter.  Its  value  was  increased  until  the  exper¬ 
imentally  measured  maximum  stress  was  obtained  (Ref.  14).  The  matrix  modulus  was 
considered  a  fixed  quantity  since  it  was  measured.  Void  properties  were  set  to  zero.  A  value 
of  Y  —  1.26(1  —  c‘)/24  was  used  in  the  M-T  models  since  this  gave  the  correct  modulus 
predictions  in  Sec.  6.3. 

Figure  11  shows  a  comparison  of  the  experimental  data  for  a  composite  containing 
clQ  =  0.30  of  untreated  beads  with  the  F-F  model  and  the  M-T  void  and  M-T  vacuole 
models.  The  parameters  used  for  these  predictions  may  be  found  in  Table  II.  It  can  be  seen 
that  the  M-T  vacuole  model  gave  the  best  prediction  for  composite  behavior.  The  M-T 
void  model  predicted  a  greater  reduction  in  modulus  and  as  a  consequence  lower  stress  at 
a  strain  >  0.25.  The  F-F  model  predicted  a  yield  point  and  a  much  lower  stress  at  a  strain 


>  0.20. 
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For  a  composite  containing  c*  =  0.40  of  untreated  beads  (Fig.  12),  tlie  F-F  model 
predicted  the  experimental  data  better.  Both  the  M-T  void  and  M-T  vacuole  models  over¬ 
predicted  the  composite  stress  after  the  yield  point.  The  behavior  of  the  composite  con¬ 
taining  4  =  0.50  of  untreated  beads  (Fig.  13)  was  more  accurately  predicted  by  the  M-T 
models  than  the  F-F  model.  It  can  be  seen  that  the  initial  modulus  predicted  by  the  M-T 
models  matched  the  experimental  data  very  well.  The  predictions  after  strain  >  0.05  cannot 
be  confirmed  because  the  specimen  ruptured  at  that  point. 

Comparison  between  the  experimental  data  for  the  composite  containing  cl0  equal  to 
0.30,  0.40  and  0.50  of  treated  beads  and  the  models  shows  that  the  M-T  vacuole  model  gave 
the  closest  predictions  (see  Figs.  14,  15  and  16).  Like  the  untreated  bead  composites  the  M- 
T  void  model  underpredicted  composite  stresses  slightly.  The  composite  containing  c'0  =  0.4 
treated  beads  did  not  have  the  large  drop  in  stress  like  its  equivalent  with  untreated  beads. 
As  seen  from  the  previous  figures,  the  M-T  models  tended  to  predict  this  type  of  behavior. 
The  F-F  model  predicted  a  yield  point  and  a  large  drop  in  stresses  afterwards.  Examination 
of  Table  II  shows  that  larger  adhesion  energies  Gc  were  needed  in  the  models  for  treated 
beads  to  account  for  the  improved  adhesion  due  to  the  silane  treatment.  These  energy 
values  provide  a  means  for  evaluating  the  adhesion  characteristics  between  the  constituents 
without  having  to  measure  it  experimentally. 

7.3  Glass  Bead/ Polyethylene  Composite 

In  Ref.  14,  the  author  fabricated  glass  bead/polyethylene  composites  with  approx¬ 
imate  inclusion  volume  fractions  of  cj,  =  0.2  and  cj,  =  0.5.  As  with  Ref.  30,  one  set  of 
composites  contained  untreated  glass  beads  while  the  other  set  used  glass  beads  treated 
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with  a  silane  coupling  agent.  The  six  composite  combinations  studied  in  Ref.  14  will  be 
used  to  demonstrate  the  performance  of  the  models.  The  matrix  modulus  was  considered 
a  fixed  quantity  for  the  M-T  models  while  the  interaction  factor  multiplier,  Ym,  had  to  be 
considered  an  adjustable  parameter  for  this  material  system.  This  will  be  discussed  shortly. 
Void  properties  were  set  to  zero.  Adhesion  energy  was  considered  an  adjustable  param¬ 
eter  for  all  models.  It  was  increased  until  experimentally  measured  maximum  stress  was 
obtained. 

Figure  17  shows  a  comparison  of  the  experimental  data  for  a  composite  containing 
c*  =  0.19  of  31  jum  untreated  beads  (designated  U2520)  with  the  F-F,  M-T  void  and  M-T 
vacuole  models.  The  F-F  model  predicted  a  yield  point  and  then  a  drop  in  stress.  The 
M-T  models  provided  better  predictions  when  Ym  =  0.6  was  used.  In  comparison  with 
Ym  =  1.26  for  glass/polyurethane  composites  (Sec.  6.3)  where  it  was  suggested  that  the 
value  of  1.26  indicated  the  presence  of  additional  reinforcement,  the  value  of  0.6  supports 
the  idea  that  there  was  an  absence  of  additional  reinforcement  due  to  the  poor  adhesion 
between  particle  and  matrix.  It  can  be  seen  that  the  M-T  vacuole  model  predicted  the 
behavior  of  the  composite  better  for  strains  >  0.025. 

The  results  for  the  composite  containing  c*  =  0.22  of  31  ^m  treated  beads  (T2520) 
shows  that  all  models  had  difficulty  with  the  nonlinearity  at  the  beginning  of  the  stress- 
strain  curve  (Fig.  18).  This  was  caused  by  the  use  of  a  linear  elastic  matrix  modulus  in  the 
models. 


The  results  for  T2520  suggest  that  the  nonlinear  effects  due  to  matrix  nonlinearity 
were  as  significant  as  the  nonlinear  effects  due  to  inclusion  debonding.  Two  factors  allowed 
this  to  happen.  First,  T2520  contained  alow  volume  fraction  of  inclusions.  This  permitted 
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matrix  properties  to  have  a  greater  influence  on  overall  behavior.  Second,  the  actual  stress- 
strain  curve  for  polyethylene  has  a  concave-down  shape  up  to  its  maximum  stress.  This 
meant  that  the  actual  modulus  was  high  initially  and  then  dropped  off  as  maximum  stress 
was  reached.  The  predicted  values  were  based  on  an  average  matrix  modulus  which  was 
defined  as  the  secant  modulus  measured  at  90%  of  maximum  stress  (Ref.  14).  Thus, 
composite  stresses  were  underpredicted  initially  because  an  average  linear  elastic  matrix 
modulus  was  used. 

The  M-T  models  for  T2520  used  a  Ym  =  2.5.  Following  the  reasoning  given  pre¬ 
viously,  this  indicates  that  the  silane  coupling  agent  was  causing  additional  reinforcement 
effects  to  be  seen.  The  difference  in  Ym  for  the  U2520  and  T2520  predictions  is  reflected 
in  the  differences  found  in  their  overall  stress-strain  behavior.  T2520  has  a  larger  ini¬ 
tial  modulus  and  a  higher  maximum  stress  than  U2520.  The  untreated  and  treated  glass 
bead/polyurethane  composites  (cj,  =  0.40)  have  similar  stress-strain  behavior  in  spite  of 
the  fact  that  one  has  a  surface  treatment.  This  explains  why  the  glass  bead/polyethylene 
composites  need  an  adjustable  Ym  while  the  glass  bead/polyurethane  composites  do  not. 

Figure  19  for  a  composite  containing  cj,  =  0.48  of  31  fx m  untreated  beads  (U2550) 
shows  that  all  models  could  reproduce  the  experimental  data.  However,  for  the  F-F  model 
to  approximate  the  experimental  behavior,  the  measured  matrix  modulus  of  187  MPa  had 
to  be  reduced  to  90  MPa.  The  justification  given  in  Ref.  14  for  this  reduction  was  that 
poor  adhesion  resulted  in  loss  of  reinforcement  and  lack  of  strain  in  the  matrix.  As  with 
the  U2520  composite,  the  M-T  models  accounted  for  this  behavior  through  modification  of 
Ym.  To  reproduce  the  experimental  data  Ym  =  0.8  was  used.  As  before,  this  value  indicates 
that  there  was  an  absence  of  additional  reinforcement  due  to  poor  adhesion. 
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A  comparison  of  the  results  for  a  composite  containing  c*0  =  0.49  of  31  /am  treated 
beads  (T2550)  in  Fig.  20  shows  that  the  M-T  models  gave  the  best  prediction  of  mechanical 
behavior.  Ym  =  2.5  was  used.  As  with  the  T2520  M-T  results,  this  value  indicated  that 
additional  reinforcement  was  being  produced  by  the  coupling  agent.  Nonlinearity  due  to 
actual  matrix  behavior  was  less  evident  for  this  composite  because  there  is  less  matrix 
volume.  This  allowed  nonlinearity  due  to  inclusion  debonding  to  dominate.  The  F-F  model 
predicted  a  yield  point  and  a  decline  in  stress  after  a  strain  of  0.075. 

Figures  21  and  22  shows  a  comparison  of  results  for  composites  containing  cla  —  0.19 
and  c*  =  0.49  of  130/^m  treated  beads  respectively.  The  c‘  —  0.19  T1020  composite  (Fig. 
21)  behaves  much  like  the  T2520  composite  (Fig.  17).  Again  the  M-T  and  F-F  models 
underpredict  the  stresses  at  strains  below  0.035.  The  reasons  given  during  the  discussion 
of  the  T2520  results  would  apply  to  the  T1020  results  as  well.  At  strains  >  0.035,  it  can 
be  seen  that  the  M-T  vacuole  results  follow  more  closely  the  experimental  results.  For  the 
cl0  —  0.49  T1050  composite  results  (Fig.  22),  the  M-T  predictions  follow  the  experimental 
data  better.  The  value  of  Ym  =  2.5  indicated  that  additional  reinforcement  was  present.  An 
examination  of  Table  III  shows  that  larger  adhesion  energies  were  required  in  the  models 
to  reflect  the  improved  adhesion  due  to  the  silane  treated  beads.  Adhesion  energies  could 
not  be  verified  because  they  were  not  measured. 

8.0  SUMMARY 

The  theoretical  framework  for  calculating  composite  modulus  using  an  improved 
Mori- Tanaka  method  has  been  presented.  A  comparison  of  2-phase  composite  results  with 
experimental  data  showed  that  the  M-T  solution  with  no  particle  interaction  effects  un- 
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derpredicted  the  measured  moduli.  The  inclusion  of  particle  interaction  effects  through  a 
correction  matrix  [Tr]  improved  predictions  but  had  some  undesirable  side  effects.  These 
side  effects  were  eliminated  by  making  the  interaction  factor  Y  a  function  of  the  inclusion 
volume  fraction. 

A  comparison  with  3-phase  modulus  data  obtained  from  the  literature  showed  there 
was  good  agreement  between  the  3-phase  M-T  model  and  the  Farber-Farris  model  for  low 
initial  inclusion  volume  fractions.  At  higher  initial  inclusion  volume  fractions,  predictions 
for  a  hypothetical  3-phase  composite  showed  that  the  3- phase  M-T  model  predicted  a  faster 
decrease  in  modulus  than  the  F-F  model.  Moduli  solutions  for  composites  containing  vac¬ 
uoles  showed  that  a  slightly  higher  modulus  could  be  expected  when  compared  to  composites 
containing  voids. 

Examination  of  the  Poisson  ratio  results  for  3-phase  composites  showed  that  a  non¬ 
monotonic  behavior  was  predicted  as  inclusion  volume  fraction  was  increased.  This  behavior 
was  restricted  to  modulus  prediction  models  which  were  derived  using  a  discrete  averaging 
process  to  account  for  the  reinforcement  of  additional  inclusions. 

Based  on  the  improved  M-T  method,  new  micromechanical  models  for  the  predic¬ 
tion  of  particulate  composite  mechanical  behavior  were  developed.  Comparisons  between 
the  M-T  void  and  M-T  vacuole  models  and  the  experimental  data  for  glass /polyurethane 
and  glass/polyethylene  composites  showed  that  the  M-T  vacuole  model  gave  the  best  re¬ 
sults  for  inclusion  volume  fractions  ranging  from  cl0  =  0.2  to  cla  =  0.4.  This  suggests  vacuole 
formation  rather  than  void  formation  is  more  representative  of  the  actual  debonding  pro¬ 
cess.  The  composites  containing  cj,  =  0.5  particles  showed  that  the  M-T  models  generally 
performed  better  than  the  F-F  model.  The  M-T  models  either  predicted  the  initial  modulus 
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or  the  stress  after  the  knee  of  the  stress-strain  curve  more  closely. 

The  introduction  of  an  interaction  factor  multiplier  Ym  in  the  M-T  models  provided 
a  means  of  evaluating  the  degree  of  particle  interaction  in  the  glass/polyethylene  composites. 
A  low  Ym  suggested  that  there  was  little  interaction  due  to  poor  adhesion  between  the 
phases.  A  high  Ym  suggested  there  was  greater  interaction  due  to  improved  adhesion. 

The  only  case  where  the  F-F  model  performed  better  than  the  M-T  models  was 
in  a  glass/polyurethane  composite  which  contained  c'0  —  0.4  of  untreated  beads.  The 
experimental  data  showed  a  yield  point  and  a  large  drop  in  stress  afterwards.  Since  the  F-F 
predictions  tend  to  have  this  shape,  the  F-F  model  gave  good  results  for  this  situation.  The 
M-T  models  tended  to  show  a  much  more  gradual  decline  in  stress.  From  the  composite 
systems  studied,  the  gradual  decline  behavior  appeared  to  be  more  common. 

The  F-F  and  M-T  models  had  trouble  predicting  the  mechanical  behavior  of  a 
glass/polyethylene  composite  containing  c'0  =  0.22  of  treated  glass  beads.  The  discrepancy 
between  theoretical  and  experimental  results  was  attributed  to  the  use  of  linear  elastic 
matrix  properties.  This  shortcoming  indicates  that  improvements  to  M-T  void  and  M-T 
vacuole  predictions  can  be  made  by  reformulation  of  the  models  to  include  nonlinear  matrix 


behavior. 
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TABLE  I 


Elastic  properties  of  constituent  phases  for  particulate  composites 


Ref. 

Matrix 

Inclusion 

G0 

(MPa) 

Gi 

(GPa) 

VQ 

Vi 

14 

polyethylene 

glass  bead 

187 

30.2 

0.34 

0.16 

25 

cobalt 

tungsten-carbide 

79000 

293 

0.31 

0.194 

29 

epoxy 

glass  bead 

1080 

30.9 

0.394 

0.23 

30 

polyurethane 

glass  bead 

1.4 

30.2 

0.499 

0.16 

35 

epoxy 

sand 

725 

29.4 

0.4 

0.25 

Go,  matrix  shear  modulus,  Gi,  inclusion  shear  modulus,  ua,  matrix  Poisson  ratio,  v,,  inclusion  Poisson 
ratio. 
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TABLE  II 


Model  input  parameters  for  glass  bead/polyurethane  composites 


Surface  Treatment 
Volume  Fraction 

Untreated 

0.3 

Untreated 

0.4 

Untreated 

0.5 

Treated 

0.3 

Treated 

0.4 

Treated 

0.5 

Avg.  rad.  (/i m) 

12.5 

12.5 

12.5 

12.5 

12.5 

12.5 

Log  std.  dev. 

0.228 

0.228 

0.228 

0.228 

0.228 

0.228 

0.30 

0.40 

0.50 

0.30 

0.40 

0.50 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Gi  (GPa) 

30.17 

30.17 

30.17 

30.17 

30.17 

30.17 

G0  (MPa) 

1.40 

1.40 

1.40 

1.40 

1.40 

1.40 

Vi 

0.16 

0.16 

0.16 

0.16 

0.16 

0.16 

Vo 

0.499 

0.499 

0.499 

0.499 

0.499 

0.499 

Pl 

yb 

0.7 

0.7 

0.7 

0.7 

0.7 

0.7 

1.26 

1.26 

1.26 

1.26 

1.26 

1.26 

GC:F-FC  (J/m2) 

8 

6 

6 

12 

8 

10 

GC:M-Td  void  (J/m2) 

12 

8 

13 

16 

14 

21 

GC:M-Te  vacuole  (J/m2) 

12 

8  . 

14 

17 

14 

24 

a  Maximum  packing  fraction  used  with  F-F  model, 
b  Interaction  factor  multiplier.  Found  by  numerical  trial  and  error, 
c  Adhesion  energy  used  with  F-F  model.  Found  by  numerical  trial  and  error, 
d  Adhesion  energy  used  with  M-T  void  model.  Found  by  numerical  trial  and  error, 

e  Adhesion  energy  used  with  M-T  vacuole  model.  Found  by  numerical  trial  and  error. 

Co,  initial  inclusion  volume  fraction,  c^,  initial  void  or  vacuole  volume  fraction,  Gi,  inclusion  shear  modulus, 
G0,  matrix  shear  modulus,  i/;,  inclusion  Poisson  ratio,  v0,  matrix  Poisson  ratio. 
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TABLE  III 


Model  input  parameters  for  glass  bead/polyethylene  composites 


Surface  Treatment 
Composite  No. 

Untreated 

U2520 

Untreated 

U2550 

Treated 

T2520 

Treated 

T2550 

Treated 

T1020 

Treated 

T1050 

Avg.  rad.  (/xm) 

15.5 

15.5 

15.5 

15.5 

65 

65 

Log  std.  dev. 

0.167 

0.167 

0.167 

0.167 

0.0374 

0.0374 

0.19 

0.48 

0.22 

0.49 

0.19 

0.49 

Co 

0.03 

0.01 

0.0 

0.0 

0.0 

0.0 

Gi  (GPa) 

30 

30 

30 

30 

30 

30 

Go : F-F  (MPa) 

187 

90 

187 

187 

187 

187 

G0:M-T  (MPa) 

187 

187 

187 

187 

187 

187 

Vi 

0.16 

0.16 

0.16 

0.16 

0.16 

0.16 

v0 

0.34 

0.34 

0.34 

0.34 

0.34 

0.34 

yl 

0.6 

0.6 

0.6 

0.6 

0.6 

0.6 

0.6 

0.8 

2.5 

2.5 

2.5 

2.5 

Gc: F-Fc  (J/m2) 

5 

2.5 

20 

6 

35 

6.0 

GC:M-Td  void  (J/m2) 

4.5 

2 

19 

12 

35 

10 

GC:M-Te  vacuole  (J/m2) 

4 

1.8 

17 

11 

33 

10 

a  Maximum  packing  fraction  used  with  F-F  model, 
b  Interaction  factor  multiplier.  Found  by  numerical  trial  and  error, 
c  Adhesion  energy  used  with  F-F  model.  Found  by  numerical  trial  and  error, 
d  Adhesion  energy  used  with  M-T  void  model.  Found  by  numerical  trial  and  error, 
e  Adhesion  energy  used  with  M-T  vacuole  model.  Found  by  numerical  trial  and  error. 

initial  inclusion  volume  fraction,  Co ,  initial  void  ox  vacuole  volume  fraction,  G , ,  inclusion  shear  modulus, 
G0:F-F,  matrix  shear  modulus  used  with  F-F  model,  G0:M-T,  matrix  shear  modulus  used  with  M-T  models, 
Vi,  inclusion  Poisson  ratio,  v0,  matrix  Poisson  ratio. 
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FIGURE  1  -  Normalized  modulus  vs.  volume  fraction  for  2-phase  glass  bead/epoxy 
composite.  Experimental  data  from  Smith  (Ref.  29). 
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FIGURE  2  -  Normalized  modulus  vs.  volume  fraction  for  2-phase  tungsten- 
carbide/cobalt  cermet.  Experimental  data  from  Ravichandran  (Ref. 
25). 
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FIGURE  3  -  Normalized  modulus  vs.  volume  fraction  for  2-phase  glass  bead/poly- 
urethane  composite.  Experimental  data  from  Yilmazer  (Ref.  30). 
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FIGURE  5  -  Composite  modulus  vs.  volume  fraction  for  3-phase  glass  bead/poly- 
urethane/void  or  vacuole  composite  using  properties  from  Ref.  30. 
cv  -  0.5  -  c\ 
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FIGURE  6  -  Composite  modulus  vs.  volume  fraction  for  3-phase  glass  bead/ 
epoxy /void  or  vacuole  composite  using  properties  from  Ref.  29. 
cv  =  0.5  -  c*. 


Inclusion  Volume  Fraction  (c>) 
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FIGURE  7  -  Composite  Poisson  ratio  vs.  volume  fraction  for  3-phase  glass  bead/ 
polyurethane/void  composite  using  properties  from  Ref.  30.  cv  = 
0.5 -cL 


Inclusion  Volume  Fraction  (c>) 
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FIGURE  9  -  Stress-strain  and  inclusion  volume  fraction  behavior  for  a  generic  glass 
bead/epoxy  composite  containing  various  levels  of  initial  inclusion  vol¬ 
ume  fractions. 
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FIGURE  10  -  Correspondance  between  inclusion  volume  fraction  and  probability  of 
survival  values  as  function  of  composite  strain. 
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FIGURE  11  -  Mechanical  behavior  predictions  for  composite  containing  c'  =  0.3  un¬ 
treated  glass  beads  in  polyurethane  (Ref.  30).  cv  =  cl0  —  cl.  Interaction 
between  bonded  particles  taken  into  account. 
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FIGURE  12  -  Mechanical  behavior  predictions  for  composite  containing  c*  =  0.4  un¬ 
treated  glass  beads  in  polyurethane  (Ref.  30).  cv  —  c!0-c*.  Interaction 
between  bonded  particles  taken  into  account. 


Strain  (mm/mm) 
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FIGURE  13  -  Mechanical  behavior  predictions  for  composite  containing  cl0  =  0.5  un¬ 
treated  glass  beads  in  polyurethane  (Ref.  30).  cv  —  c'0  —  c*.  Interaction 
between  bonded  particles  taken  into  account. 


Strain  (mm/mm) 
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FIGURE  14  -  Mechanical  behavior  predictions  for  composite  containing  c*  =  0.3 
treated  glass  beads  in  polyurethane  (Ref.  30).  cv  =  c'-c*.  Interaction 
between  bonded  particles  taken  into  account. 


Strain  (mm/mm) 
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FIGURE  15  -  Mechanical  behavior  predictions  for  composite  containing  c'0  —  0.4 
treated  glass  beads  in  polyurethane  (Ref.  30).  cv  —  cl0  —  cl.  Interaction 
between  bonded  particles  taken  into  account. 
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FIGURE  16  -  Mechanical  behavior  predictions  for  composite  containing  c\  -  0.5 
treated  glass  beads  in  polyurethane  (Ref.  30).  cv  —  cl0  -  cl .  Interaction 
between  bonded  particles  taken  into  account. 


Strain  (mm/mm) 
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FIGURE  17  -  Mechanical  behavior  predictions  for  composite  containing  cl0  —  0.19 
31  jam  untreated  glass  beads  in  polyethylene  (Ref.  14).  cv  =  c‘0  -  c*. 
Interaction  between  bonded  particles  taken  into  account. 
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FIGURE  18  -  Mechanical  behavior  predictions  for  composite  containing  c*  =  0.22 
31  fim  treated  glass  beads  in  polyethylene  (Ref.  14).  cv  =  c‘0  -  c‘. 
Interaction  between  bonded  particles  taken  into  account. 
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FIGURE  19  -  Mechanical  behavior  predictions  for  composite  containing  c’  =  0.48 
31  n m  untreated  glass  beads  in  polyethylene  (Ref.  14).  cv  —  c{0  —  c*. 
Interaction  between  bonded  particles  taken  into  account. 
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Mechanical  behavior  predictions  for  composite  containing  c’  =  0.49 
31  nm  treated  glass  beads  in  polyethylene  (Ref.  14).  cv  =  c*  -  c\ 
Interaction  between  bonded  particles  taken  into  account. 


Strain  (mm/mm) 
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FIGURE  21  -  Mechanical  behavior  predictions  for  composite  containing  cl0  =  0.19 
130  yum  treated  glass  beads  in  polyethylene  (Ref.  14).  cv  =  c%0  —  c*. 
Interaction  between  bonded  particles  taken  into  account. 
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FIGURE  22  -  Mechanical  behavior  predictions  for  composite  containing  cl0  =  0.49 
130  pm  treated  glass  beads  in  polyethylene  (Ref.  14).  cv  =  clQ  —  c*. 
Interaction  between  bonded  particles  taken  into  account. 


Strain  (mm/mm) 
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Mathematica  Listing  for  3-Phase  Modulus  Model 


Program  Listing 

(*  Three-phase  Mori-Tanaka  model  with  particle  interaction  effects  \ 
and  void  or  vacuole  formulation  *) 

(*  cx  -  inclusion  volume  fraction,  \ 

cy  -  void  or  vacuole  volume  fraction,  \ 
jucor  -  apply  ju  correction  (l=yes) ,  \ 
voidmag  -  apply  correction  to  voids  (l=yes)  *) 

cx=0.30119;  cy=0.5-cx;  jucor=l;  voidmag=0 

(*  particle  interaction  function  \ 

ygr  -  inclusion  interaction  function,  \ 
ygs  -  void  or  vacuole  interaction  function  *) 

ygr=1.26*(l-cx)/24;  ygs=(l)/24 

(*  material  constants  for  comparison  and  inclusion  materials  *) 

(*  phases  in  terms  of  tensile  modulus  and  Poisson’s  ratio  *) 

(*  em  -  matrix  tensile  modulus,  \ 
v  -  matrix  Poisson  ratio,  \ 
emr  -  inclusion  tensile  modulus,  \ 
vr  -  inclusion  Poisson  ratio,  \ 
ems  -  void  or  vacuole  tensile  modulus,  \ 
vs  -  void  or  vacuole  Poisson  ratio,  \ 
k  -  matrix  bulk  modulus,  \ 
u  -  matrix  shear  modulus,  \ 

es_ij  -  vacuole  tensile  modulus  in  ij -direction,  \ 
vs_ij  -  vacuole  Poisson  ratio  in  ij-direction  *) 

(*  Smith  glass  bead/epoxy  properties  *) 

(*  em=3.01  ;  v=0.394;  \ 
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emr='76  ;  vr=0.23;  \ 

ems=0  ;  vs=0.23;  \ 

k  =  em/(3*(l-2*v))  ;  \ 

u  =  em/(2*(l+v)) ;  \ 

kr  =  emr/(3*(l-2*vr)) ;  \ 

ur  =  emr/(2*(i+vr)) ;  \ 

ks  =  ems/(3*(l-2*vs)) ;  \ 

us  =  ems/(2*(l+vs)) ;  \ 


esll=0.0*enis;  es22=es33=i . 0*ems ;  \ 
vsi2=vsl3=0;  \ 
vs21~vs31=vs23=vs32=vs  *) 

(*  Ravichandran  #625  *) 


(*  em=207  ;  v=0.31;  \ 

emr=700  ;  vr=0.194;  \ 

eros=0  ;  vs=0.194;  \ 

k  =  em/(3*(l-2*v)) ;  \ 

u  =  em/(2*(l+v)) ;  \ 

kr  =  emr/(3*(i-2*vr)) ;  \ 

ur  =  emr/(2*(l+vr) ) ;  \ 

ks  =  ems/(3*(i-2*vs) ) ;  \ 

us  =  ems/(2*(l+vs) ) ;  \ 


esll=l .0*ems ;  es22=es33=l . 0*ems ;  \ 
vsl2=vsl3=0;  \ 
vs21=vs31=vs23=vs32=vs  *) 


(*  Yilmazer  #351  *) 
em=4.2  ;  v=0.499;  \ 

emr=70000  ;  vr=0.16;  \ 

ems=0000  ;  vs=0.16;  \ 

k  =  em/(3*(l-2*v)) ;  \ 

u  =  em/(2*(l+v) ) ;  \ 

kr  =  emr/(3*(l-2*vr)) ;  \ 

ur  =  emr/(2*(l+vr)) ;  \ 

ks  =  ems/(3*(l-2*vs)) ;  \ 

us  =  ems/(2*(l+vs)) ;  \ 


esll^O.O+ems;  es22=es33=l .0*ems ;  \ 
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vsl2=vsl3=0;  \ 
vs21=vs31=vs23=vs32=vs 

(*  interaction  parameters  for  inclusions  *) 
alphar  =  2*(5*v-l)+10*(i-v)*(k/(kr-k)-u/(ur-u)) 
betar  =  2*(4-5*v)+i5*(i-v)*(u/(ur-u)) 

zet  air  =  12*v*  ( 13-  14*v)  -  (96*alphar/  (3*alphar+2*betar ) )  *  (  i-2*v)  *  ( 1+v) 
zeta2r  =  6*(25-34*v+22*v~2)-(36*alphar/ (3*alphar+2*betar))*(l-2*v)*(l+v) 
magr  =  ((5/4)/betar~2)*ygr 

(*  interaction  parameters  for  voids  *) 
alphas  =  2*(5*v-l)+10*(l-v)*(k/(ks-k)-u/(us-u)) 
betas  =  2*(4-5*v)+15*(l-v)*(u/(us-u)) 

zetals  =  12*v* ( 13-14*v) - (96*alphas/ (3*alphas+2*betas) ) * ( l-2*v) * ( 1+v) 
zeta2s  =  6*(25-34*v+22*v~2)-(36*alphas/(3*alphas+2*betas))*(l-2*v)*(l+v) 
mags  =  ( (5/4)/betas~2)*ygs 

(*  Eshelby’s  transformation  factors  for  spherical  particles  *) 
si  =  5*v-l 
s2  =  4-5*v 
sdet  =  15*(l-v) 

C*  comparison  or  matrix  material  stiffness  matrix  *) 
ciill  =  k+(4/3)*u;  cll22  =  k-(2/3)*u;  cll33  =  cll22 
c2211  =  C1122  ;  c2222  =  cllll  ;  c2233  =  cll22 

C3311  =  cll22  ;  c3322  =  cll22  ;  c3333  =  cllll 

cO  =  ■C{cllll,cll22,cll33>,{c2211,c2222,c2233>,{c3311,c3322,c3333» 

(*  inclusion  material  stiffness  matrix  *) 

crllll  =  kr+(4/3)*ur;  crll22  *  kr-(2/3)*ur;  crll33  =  crll22 
cr2211  =  crll22  ;  cr2222  =  crllll  ;  cr2233  =  crll22 

cr3311  =  crll22  ;  cr3322  =  crll22  ;  cr3333  =  crllll 

cr  =  {-Ccrllll  ,crll22  ,crll33>  ,-Ccr2211,cr2222,cr2233}J'Ccr3311,cr3322,cr3333» 


(*  void  or  vacuole  material  stiffness  matrix  *) 
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csdet  =  1-vs12*vs21-vs13*vs31-vs23*vs32-vs12*vs23*vs31-vs13*vs21*vs32 

csllll  =  esll(l-vs23*vs32) 

csll22  “  esll(vs21+vs23*vs31) 

csll33  =  esll(vs31+vs2i*vs32) 

cs2211  =  es22(vsl2+vsl3*vs32) 

cs2222  =  es22(l-vsl3*vs31) 

cs2233  =  es22(vs32+vsl2*vs31) 

cs3311  =  es33(vsl3+vsl2*vs23) 

cs3322  =  es33(vs23+vsl3*vs21) 

cs3333  =  es33(l-vsl2*vs21) 

cs  =  (l/csdet)*{{csllll ,csll22,csll33},  \ 

{cs2211,cs2222,cs2233},  \ 

{cs33il ,  cs3322 ,  cs3333» 

(*  Eshelby  transformation  tensor  for  spherical  particles*) 

sllil  =  si+2*s2 ;  sll22  =  si  ;  sll33  =  sll22 

s221:L  =  sll22  ;  s2222  =  sllil  ;  s2233  =  sll22 

s331i  =  S1122  ;  s3322  =  s!122  ;  s3333  =  sllil 

s  =  (l/sdet)*{-[sllll,sll22,sll33>,  \ 

-Cs2211,s2222,s2233>,  \ 

{s3311 ,  s3322 ,  s3333» 


(*  correction  factor  matrix  for  inclusions  *) 
cwrll  =  zetalr+2*zeta2r  ;  cwrl2  =  zetalr  ;  cwrl3  =  cwrl2 

cwr21  =  cwrl2  ;  cwr22  =  cwrli  ;  cwr23  =  cwrl2 

cwr3:L  =  cwrl2  ;  cwr32  =  cwrl2  ;  cwr33  =  cwrll 

cwr  ~  {{cwrll ,cwrl2,cwrl3>,  \ 

{cwr21,cwr22,cwr23>,  \ 

{cwr3 1 , cwr32 , cwr33L} 


C*  correction  factor  matrix  for 


cwsll  =  zetals+2*zeta2s 
cws21  =  cwsl2 
cws3:l  =  cwsl2 


cwsl2 

cws22 

cws32 


voids  *) 
=  zetals 
=  cwsll 
=  cwsl2 


;  cwsl3  =  cwsl2 
;  cws23  =  cwsl2 
;  cws33  -  cwsll 
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cws  =  {{cwsll ,cwsl2,cwsl3>,  \ 

{cws2i ,cws22,cws23},  \ 

{cws3 1 ,  cws32 ,  cws33}} 

(*  Identity  matrix  *) 
i  =  IdentityMatrix[3] 

(*  A-matrices  *) 

ar  =  Inverse [cr-cO] .cO 

as  =  Inverse [cs-cO] . cO 

(*  complete  correction  matrix  for  inclusions  and  voids  *) 

If [ jucor>0 , lambdar=i+magr*cx*cwr , lambdar=i] 

If [voidmag>0 , lambdas=i+mags*cy*cws , lambdas=i] 

(*  calculation  of  average  stiffness  coefficients  *) 
fl  =  Inverse [s+ar] 
f2  =  Inverse [s+as] 

f 3  =  Inverse [cx* (i-lambdar-s)+s+ar+cy* (i-lambdas-s) .f 2. (s+ar)] 
f4  =  Inverse [cy*(i-lambdas-s)+s+as+cx*(i-lambdar-s) .f 1 . (s+as)] 
fcavg  =  cO. (i+cx*lambdar .f3+cy*lambdas .f4) 

f emav  =  fcavg [[1,1]] - (2*f cavg [[1,2]] *f cavg [[2,1]])/ (fcavg [[2,2]] +f  cavg [[2,3]]) 
fvav  =  f  cavg [ [2 , 1] ] / (f  cavg [ [2 , 2] ] +f  cavg [ [2 , 3]  ]  ) 

Print["cx=  ",cx] ;  Print["cy=  ",cy];  \ 

Print["ju  corrc=  " , jucor] ;Print["void  mag=  ".voidmag];  \ 

Print ["f  eavg=  " ,f emav] ;  Print ["f  vavg=  M,fvav] 

(*  use  for  generation  of  modulus  and  Poisson  ratio  vs  inclusion  Vf  *) 

(*  fe  =  N [Table [f emav, {cx, 0,0. 5, 0.01}] ,4]  \ 
fv  =  N[Table[fvav,{cx,0,0.5,0.01}] ,4]  *) 
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APPENDIX  B 


FORTRAN  Listing  for  Micromechanical  Model 


Program  Listing 

C====  main  program 

C  calculates  particle  size  histogram  with  corresponding  filler 
C  volume  fraction,  uses  Z-decrements  for  particle  size 
C  determination 

C  user  enters  the  following  information: 

C  number  of  points  desired  in  stress-strain  curve 

C  number  of  particle  distributions 

C  avg  radius  and  std  dev  of  each  distribution 

C  volume  fraction  of  filler  and  voids  of  each  distribution 

C  matrix  and  filler  shear  modulus,  poisson  ratio 

G  void  shear  and  bulk  moduli 

C  sample  volume,  fraction  debond,  w-type,  m-type 

C  adhesion  energy  and  applied  pressure 

C  information  may  be  entered  using  keyboard  or  by  input  data  file. 

C  three  options  for  printing  out  intermediate  results  are  available: 
C  if  values  for  no  pts  desired  in  stress-strain  curve  and  number 

C  of  particle  distributions  are  negative,  data  files  GAUSS, 

G  HISTO,  DEBUG  and  STRESS  are  written. 

C  if  value  for  no  pts  desired  in  stress-strain  curve  is  negative 

C  and  value  for  number  of  particle  distributions  is  positive, 

C  data  files  HISTO,  DEBUG  and  STRESS  are  written. 

C  if  input  data  was  entered  using  data  file,  the  data  file  STRESS 
C  will  be  renamed  to  the  input  data  file’s  name. 

C 

C  to  compile  and  link:  fl  p.for  /link  graphics .lib .  the  files 

C  MSGRAPH.FOR  and  GRFDEF.FOR  should  be  in  the  same  directory  unless 

C  a  temporary  variable  has  been  set  up  to  point  to  the  location  of 

C  include  files,  these  files  contain  graphics  routines  necessary  to 
C  plot  stress-strain  curve  on  screen. 

C 

C  the  indexing  of  the  various  parameters  is  organized  as  follows : 
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C  strain  -  current  critical  strain 

C  stress  -  current  stress  calc  using  previous  E,G  and  crit. strain 

C  Pr_surv  -  current  no.  of  particles  remaining 

C  E,G  -  moduli  at  current  Pr.surv 

C  Vf  ,Vv  -  current  filler  and  void  volume  fractions 

C  dG/dc,  dK/dc,  dA/dc  -  current  differential  quantities 

C  the  total  number  of  points  is  NDIST*NPTS+1  where  the  additional 

C  point  is  for  zero  strain  and  stress .  the  first  group  of  debonded 
C  particles  begins  at  ICNT=2. 

C 

C  added  routine  to  output  SQRT(r*dE/dc) .  dc  based  on  total  volume 

C  instead  of  Vf+Vm.  trapped  zero  in  SQRT  calc  of  crit.  strain. 

C 

C  implementation  of  Mori-Tanaka  solution  extended  for  3-phase  and 
C  particle  interaction,  constituent  material  properties 

C  designated  as  follows:  1-inclusion, 2-void  or  vacuole, 3-matrix. 

C  fraction  debond  (FDBND)  for  orthotropic  properties  in  loading 
C  direction,  multiplier  for  rad.  dist.  func.  (YMULT) ,  w-type 
C  designates  use  inclusion  or  void  properties  in 

C  calc  of  Wv  matrix  (0=void,  l=inclusion) ,  m-type  determines  type,  of 
C  particle  interaction  used  (O=none,i=inclusion,  2=inclusion  and 
C  void  or  vacuole) ,  v-type  determines  isotropic  or  orthotropic  matl  , 
C  (O=orthotropic,l=isotropic) 

C 

C  last  revision:  04  MAR  1995 

C 

C  set  PTMX  =  NTDIS*GSMX 

C 

REAL  LOGSTD ,  NPARTL ,  NTJMPAR ,  NETVF ,  NETVV 
REAL  IDENT , K , KCMP , MAG 
INTEGER  GSMX,PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750,NTDIS  »  3) 

COMMON  /GAUS/  Z (GSMX) .RADIUS (NTDIS, GSMX) ,PR0B(NTDIS, GSMX) 

COMMON  /DEBUG/  NUMPAR(NTDIS , GSMX) ,VOLPAR(NTDIS , GSMX) , 

*  NETVF (PTMX) .NETVV(PTMX) .DADC(PTMX) .NPARTL (NTDIS) 

COMMON  /DIST/  RADAVG(NTDIS) .LOGSTD (NTDIS) .VLFRFO (NTDIS) , 

*  VLFRVO (NTDIS) 
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COMMON  /MATRA/  BETA(2) ,WI(3,3) ,WV(3,3) ,IDENT(3,3) 

COMMON  /MATRB/  S(3,3) ,CA(3,3) ,CB(3,3) ,CE(3,3) ,CF(3,3) 

COMMON  /PROPA/  K(3) ,G(3) ,E(3) ,P0IS(3) ,CI(3,3) ,CV(3,3) ,C0(3,3) 
COMMON  /PROPB/  Cll(PTMX) ,C12(PTMX) ,C21(PTMX) ,C22(PTMX) ,C23(PTMX) , 

*  ECMP (PTMX) , P  0 ISC (PTMX ) 

COMMON  /RESULT/  CRTSTN(PTMX) .STRESS (PTMX) ,DILAT(PTMX) , 

*  PRBSRV(PTMX) ,SORRAD(PTMX) .SORPAR(PTMX) .SORVLP(PTMX) , 

*  IPDIST(PTMX) 

CHARACTER  FILNM*8 

C 

C==  initialize  variables  and  arrays  by  BLOCK  DATA  INIT 
C 

CALL  INPUT (NDIST , NTOT , VOLSMP , FDBND , YMULT , IKIND , IMORI , IPO IS , GAMM , 

*  PRESS, FILNM.IWRT) 

C 

IABORT  =0.0 

CALL  STRSTN (NDIST , NTOT , NPTS .VOLSMP , FDBND , YMULT , IKIND , IMORI , IPOIS , 

*  GAMM, PRESS, DILATO.IWRT, IABORT) 

C 

C==  write  out  results  and  debug  data 

CALL  STRWRT (NDIST , NPTS , VOLSMP , GAMM , FDBND , YMULT , IKIND , IMORI , IPO IS , 

*  PRESS, DILATO.FILNM, IABORT) 

IF  (ABS(IWRT) .GE.l)  CALL  DBGWRT(NDIST, NPTS, IABORT) 

IF  (ABS(IWRT) .GE.l)  CALL  DBGRAT(NDIST, NPTS, IABORT) 

C 

CALL  CRVPLT (NDIST, NPTS, IABORT) 

C 

END 

C 

C 

SUBROUTINE  INPUT (NDIST , NTOT .VOLSMP , FDBND , YMULT , IKIND , IMORI , IPOIS , 

*  GAMM, PRESS, FILNM.IWRT) 

C====  reads  in  problem  input  either  by  file  or  keyboard,  if  data  entered 
C  through  a  file,  user  inputs  name  only,  a  file  extension  of  DAT  is 

C  assumed,  the  first  line  in  the  input  file  is  used  for  a  user 

C  heading  and  is  not  read  in,  constituent  material  properties 

C  designated  as  follows:  1-inclusion, 2-void  or  vacuole  ,3-matrix 
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C  set  PTMX  =  NTDIS*GSMX 
C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
REAL  IDENT , K , KCKP , MAG 
INTEGER  GSMX.PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750, NTDIS  =  3) 

COMMON  /DIST/  RADAVG(NTDIS) , LOGSTD (NTDIS) ,VLFRFO(NTDIS) , 

*  VLFRVO (NTDIS) 

COMMON  /PROPA/  K(3) ,G(3) ,E(3) ,P0IS(3) ,CI(3,3) ,CV(3,3) ,C0(3,3) 
CHARACTER  ANS*1 ,FILNM*8 
C 

WRITE  (6,’(/,A)’)  5  Read  data  from  file?  (Y/N) 5 
READ  (5, ’ (Al) ’ )  ANS 
C 

IF  (ANS.EQ. 5Y’ )  THEN 

WRITE  (6, 5 (A)’)  5  File  to  read?  (.INP  will  be  appended) 5 
READ  (5, * (A8) ’ )  FILNM 

OPEN  (UNIT=7,FILE=FILNM//> .INP’ , FORM=  *  FORMATTED ’ ,STATUS=’QLD’ ) 
READ  (7,  *) 

READ  (7,*)  NTOT 
READ  (7,*)  NDIST 
DO  10  I  =  l.ABS(NDIST) 

READ  (7,*)  RADAVG(I) 

READ  (7,*)  LOGSTD(I) 

READ  (7,*)  VLFRFO (I), VLFRVO (I) 

10  CONTINUE 

READ  (7,*)  VOLSMP 

READ  (7,*)  FDBND , YMULT , IKIND , IMORI , IPO IS 
READ  (7,*)  G(3),G(1) 

READ  (7,*)  P0IS(3) ,P0IS(1) 

READ  (7, *)  G(2) ,K(2) 

READ  (7,*)  GAMM, PRESS 
CLOSE  (7) 

ELSE 

WRITE  (6,’(/,A,I3,A)’) 

*  ’  no.  pts  desired  in  stress-strain  curve  (‘C’.GSMX,’)’ 
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READ  (5,*)  NTOT 

WRITE  (6, 5 (A, II, A) ’)  ’  no.  of  particle  distributions  (<=’ , 

*  NTDIS ,  ’ )  ’ 

READ  (5,*)  NDIST 

DO  20  I  -  i,ABS (NDIST) 

WRITE  (6, ’ (A, II, A) ’)  ’  for  distribution  no.  ’ ,1, 

*  ’  mean  radius  (micron) ’ 

READ  (5,*)  RADAVG(I) 

WRITE  (6,’ (A)5)  ’  log  normal  radius  std  dev’ 

READ  (5,*)  LOGSTD(I) 

WRITE  (6, ’(A)’)  ’  initial  volume  fraction  filler  and  void’ 
READ  (5,*)  VLFRFO(I) ,VLFRV0(I) 

20  CONTINUE 

WRITE  (6, ’(A)’)  ’  sample  volume  (mm3)’ 

READ  (5,*)  VOLSMP 

WRITE  (6, ’(A)’)  ’  dbnd  frac.rad  dist  mult,w-type,m-type,v-type’ 

READ  (5,*)  FDBND , YMULT , IKIND , IMORI , IPOIS 

WRITE  (6, ’(A)’)  ’  matrix  and  filler  shear  modulus  (Pa)’ 

READ  (5,*)  G(3),G(i) 

WRITE  (6, ’(A)’)  ’  matrix  and  filler  Poisson  ratio’ 

READ  (5,*)  P0IS(3) ,P0IS(1) 

WRITE  (6 , ’ (A) ’ )  ’  void  shear  and  bulk  modulus  (Pa) ’ 

READ  (5,*)  G(2),K(2) 

WRITE  (6, ’(A)’)  ’  Gc  (J/m2)  and  applied  pressure  (Pa)’ 

READ  (5,*)  GAMM, PRESS 
FILNM  =  ’DEFAULT’ 

END  IF 
C 

C==  set  write  file  flag,  0=STRWRT,  1=STRWRT , DBGWRT , HSTWRT ,  2=all 
IWRT  =  0 

IF  (NTOT. LT.O. AND. NDIST. LT.O)  IWRT  =  2 
IF  (NTOT. LT.O. AND. NDIST. GT.O)  IWRT  =  1 
NDIST  «  ABS (NDIST) 

NTOT  =  ABS (NTOT) 

C 

RETURN 

END 
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C 

C 


SUBROUTINE  STRSTN (NDIST , NTOT , NPTS , VOLSMP , FDBND , YMULT , IKIND , IMORI , 

*  IPOIS , GAMM, PRESS ,DILATO , IWRT , IABORT) 

C====  main  subroutine  which  organizes  particle  size  distribution, 

C  composite  property,  critical  strain  and  stress  and  dilation 
C  calculation  modules . 

C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
INTEGER  GSMX , PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750,NTDIS  =  3) 

COMMON  /DEBUG/  NUMPARCNTDIS , GSMX) ,VOLPAR(NTDIS , GSMX) , 

*  NETVF (PTMX) .NETVV(PTMX) ,DADC(PTMX) , NPARTL (NTDIS) 

C 

C==  initialize  abort  flag 
IABORT  =  0 
C 

C==  create  gaussian  distribution  of  particle  size 

WRITE  (6,5(/,A)5)  5  Generating  particle  distribution* 

CALL  GAUSS (NDIST , NTOT , NPTS , IABORT) 

C==  write  out  gaussian  cumulative  data 

IF  (ABS(IWRT) .GE.2)  CALL  GAUWRT(NDIST, NPTS, IABORT) 

C 

C==  find  size  and  number  of  particles  to  debond 

WRITE  (6,’(/,A)5)  5  Finding  particle  size  and  number5 
CALL  PARTSZ (NDIST , NPTS , VOLSMP , IABORT) 

C==  write  out  particle  size  and  number  histogram 

IF  (ABS(IWRT) .GE.l)  CALL  HSTWRT (NDIST, NPTS, IABORT) 

C 

WRITE  (6,’(/,A)5)  5  Sorting  particle  distributions 5 
CALL  SORTER (NDIST, NPTS, IABORT) 

WRITE  (6, 5 (A)5)  5  Calculating  vol  fractions  and  dA/dc5 
CALL  VOLFRC (NDIST , NPTS , VOLSMP , IABORT) 

C 

WRITE  (6,’(/,A)5)  5  Generating  stress-strain  curve5 
C==  calculate  initial  composite  properties 
ICNT  =  1 
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CONCI  ■  NETVF(ICNT) 

CONCV  =  NETVV (ICNT) 

CALL  MTPEP (CONCI , CONCV , ICNT , FDBND , YMULT , IKIND , IMORI , IPOIS , I ABORT) 
IF  (IABORT.EQ.O)  WRITE  (6, ’ (A, IX, 13, A, 13, A) * ) 

*  ’  Calculating  point:  5 ,ICNT, ’/’ ,NDIST*NPTS+1 , ’  max’ 

C 

C==  main  routine  for  debonding  and  stress-strain  calculation. 

C  offset  pointer  ICNT  by  1  to  make  room  for  undebonded  state. 

C 

DO  10  ICNT  =  2 ,NDIST*NPTS+1 
CONCI  =  NETVF(ICNT) 

CONCV  =  NETVV (ICNT) 

CALL  MTPRP (CONCI , CONCV , ICNT , FDBND , YMULT , IKIND , IMORI , IPOIS , 

*  IABORT) 

IF  (IABORT.EQ.O)  WRITE  (6, ’ (A, IX, 13, A, 13, A) ’ ) 

*  ’+Calculating  point:  ’.ICNT, 5 / ’ ,NDIST*NPTS+1 , *  max’ 

CALL  CRIT ( ICNT , VOLSMP , GAMM , PRESS , IABORT) 

CALL  CALVAL (ICNT , PRESS , DILATO , IABORT) 

10  CONTINUE 
C 

RETURN 

END 

C 

C 

SUBROUTINE  GAUSS (NDIST , NTOT , NPTS , IABORT) 

C====  Program  calculates  the  cumulative  area  underneath  the 

C  gaussian  curve  between  the  limits  +/-  (IEND/FACT)  in  increments 

C  of  IDELT/FACT.  NTOT  is  used  to  calculate  an  appropriate  IDELT. 

C  since  IDELT  is  rounded  down,  the  exact  number  of  points  may  be 

C  greater,  this  is  reflected  in  NPTS. 

C  Particle  radii  converted  from  microns  to  millimeters. 

C  An  IEND  of  3301  gives  a  cumulative  distribution  which  starts 

C  at  0.0005  and  ends  at  0.9995.  This  avoids  having  extremely  large 

C  particles  when  the  log  standard  deviation  is  large. 

C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
INTEGER  GSMX.PTMX 
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PARAMETER  (GSMX  =  250.PTMX  =  750,NTDIS  =  3) 

PARAMETER  (ISTART  =  2.IEND  =  3301, FACT  =  1000,BEGNPT  =  0) 
EXTERNAL  FUNC 

COMMON  /GAUS/  Z(GSMX) , RADIUS (NTDIS, GSMX) ,PROB(NTDIS, GSMX) 

COMMON  /DIST/  RADAVG(NTDIS) .LOGSTD (NTDIS) ,VLFRFO(NTDIS) , 

*  VLFRVO (NTDIS) 

C 

IMAX  =  GSMX/2 

IDELT  =  2*INT( (IEND-ISTART) /NTOT) 

NPTS  ■  2* (INT( (IEND-ISTART)/IDELT)+1) 

IF  (NPTS. GT. IMAX)  THEN 

WRITE  (6, 5 (A) 5 )  5  Too  many  points:  SBR  GAUSS.5 

WRITE  (6, 5 (A, 14, A)5)  5  Over  max  dim  by  5 .NPTS-GSMX, 5  points.5 

I ABORT  =  1 

RETURN 

ELSE 

ENDIF 

C 

DO  20  J  =  l.NDIST 

IPTS  =  0 

DO  10  I  =  ISTART, IEND, IDELT 
IPTS  =  IPTS+1 
ENDPT  =  (FL0AT(I)-1)/FACT 
Z(NPTS/2+IPTS)  =  ENDPT 
Z(NPTS/2-IPTS+l)  =  -ENDPT 

C==  calculate  upper  portion  of  probability  curve 

RADTMP  =  10**(ALOG10(RADAVG(J) )+ENDPT*LOGSTD( J) ) 

RADIUS (J,NPTS/2+IPTS)  =  RADTMP/1000 
CALL  QSIMP (FUNC , BEGNPT , ENDPT , SURF) 

PROB(J ,NPTS/2+IPTS)  =  0.5+SURF 
C==  calculate  lower  portion  of  probability  curve 

RADTMP  =  10**(ALOG10(RADAVG(J) )-ENDPT*LOGSTD( J) ) 

RADIUS (J,NPTS/2-IPTS+l)  =  RADTMP/1000 
PROB(J .NPTS/2-IPTS+1)  =  0.5-SURF 
10  CONTINUE 
20  CONTINUE 


C 
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RETURN 

END 


C 

FUNCTION  FUNC(X) 

C====  function  used  for  gaussian  curve,  called  from  sbr  GAUSS,  sbr 
C  QSIMP  and  sbr  TRAPZD . 

PI  =  3.141592654 

FUNC  =  ( 1 . 0/SQRT (2 . 0*PI) ) *EXP (-X**2/2 . 0) 

RETURN 

END 

C 

C 

SUBROUTINE  PARTSZ (NDIST , NPTS , VOLSMP , I ABORT) 

C====  sbr  finds  the  total  particle  volume  on  a  per  particle  basis. 

C  from  this  the  number  of  particles  present  in  the  composite  is 
C  calculated  knowing  the  initial  volume  the  particles  occupy. 

C  the  incremental  probability  of  the  largest  particles  is 

C  calculated  using  a  fraction  (PFRAC)  of  the  previous  probability 

C  increment  so  that  there  is  a  smooth  transition  from  largest  to 
C  smaller  particle  sizes  in  terms  of  number. 

C 

C  set  PTMX  =  GSMX*NTDIS 
C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
INTEGER  GSMX, PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750,NTDIS  =  3) 

PARAMETER  (PI  =  3 . 1415927 ,PFRAC=0 .75) 

COMMON  /GAUS/  Z(GSMX) , RADIUS (NTDIS, GSMX) ,PROB(NTDIS, GSMX) 

COMMON  /DEBUG/  NUMPAR(NTDIS , GSMX) ,VOLPAR(NTDIS , GSMX) , 

*  NETVF (PTMX) , NETVV (PTMX) .DADC(PTMX) , NPARTL (NTDIS) 

COMMON  /DIST/  RADAVG(NTDIS) .LOGSTD (NTDIS) .VLFRFO (NTDIS) , 

*  VLFRVO (NTDIS) 

C 

IF  (IABORT.EQ . 1)  RETURN 
C 

C==  find  total  number  of  particles  in  given  filler  volume 
DO  20  IDIST  =  1, NDIST 
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VOLTOT  =  0 

C  find  total  volume  on  a  per  particle  basis 
DO  10  IPTS  =  l.NPTS 

IF  (IPTS.EQ.NPTS)  THEN 

VOLPAR (IDIST, IPTS)  =  PFRAC* 

*  CPROB (IDIST , IPTS) -PROB (IDIST , IPTS- 1 ) ) 

*  *(4.0/3. 0) *PI*RADIUS ( IDIST , IPTS) **3 
ELSE 

VOLPAR (IDIST, IPTS)  = 

*  (PROB (IDIST , IPTS+i) -PROB (IDIST , IPTS) ) * 

*  (4 . 0/3 . 0) *PI*RADIUS (IDIST , IPTS) **3 
END  IF 

VOLTOT  =  VOLTOT+VOLPAR(IDIST,IPTS) 

10  CONTINUE 

C  find  total  number  of  particles 

NPARTL (IDIST)  =  VLFRFO (IDIST) *VOLSMP/VOLTOT 
20  CONTINUE 

C 

C  calculate  volume  taken  up  by  particles  of  radius  r 
DO  30  IDIST  =  l.NDIST 
DO  40  IPTS  =  l.NPTS 

IF  (IPTS.EQ.NPTS)  THEN 

NUMPAR( IDIST, IPTS)  =  NPARTL (IDIST)* 

*  PFRAC* (PROB (IDIST , IPTS) -PROB ( IDIST , IPTS-1) ) 
IF  (NUMP AR(IDIST, IPTS) .LT. 1.0)  IFLAG  =  1 
VOLPAR(IDIST.IPTS)  =  NUMPAR(IDIST.IPTS)* 

*  (4.0/3. 0) *PI*RADIUS (IDIST , IPTS) **3 
ELSE 

NUMPAR (IDIST, IPTS)  =  NPARTL (IDIST)* 

*  (PROB (IDIST , IPTS+l) -PROB (IDIST , IPTS) ) 

IF  (NUMPAR(IDIST, IPTS) .LT. 1.0)  IFLAG  =  1 
VOLPAR (IDIST, IPTS)  =  NUMPAR (IDIST, IPTS)* 

*  (4 . 0/3 . 0) *P I* RADIUS (IDIST , IPTS) **3 
END  IF 

C 

IF  ( IFLAG. EQ.l)  THEN 
WRITE(6,5000) 
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*  IDIST,IPTS,RADIUS(IDIST,IPTS) ,NUMPAR(IDIST,IPTS) 
IFLAG  =  0 

ELSE 

ENDIF 

C 

40  CONTINUE 

30  CONTINUE 
C 

C  RETURN 

5000  FORMAT ( ’  Error  SBR  PARTSZ:  IDIST=’,I1,’  IPTS=\I3,’  RAD=’,E11.6, 

*  »  NUMPAR=’ ,E11.6) 

END 

C 

C 

SUBROUTINE  SORTER(NDIST , NPTS , IABORT) 

C=t=!==  loads  radius,  number  of  particles  and  total  volume  of  particles 
C  of  radius  r  from  each  distribution  in  a  master  array  to  sort. 

C  after  sorting  radii  in  ascending  order,  arrays  are  flipped 

C  according  to  radius  to  give  descending  order. 

C 

C  set  PTMX  =  NTDIS*GSMX 

C 

REAL  L0GSTD , NPARTL , NUMPAR , NETVF , NETVV 
INTEGER  GSMX , PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750,NTDIS  =  3) 

COMMON  /GAUS/  Z(GSMX) , RADIUS (NTDIS, GSMX) ,PROB(NTDIS, GSMX) 

COMMON  /DEBUG/  NUMPAR (NTDIS, GSMX) .VOLPAR (NTDIS, GSMX) , 

*  NETVF (PTMX) .NETVV(PTMX) ,DADC(PTMX) .NPARTL (NTDIS) 

COMMON  /RESULT/  CRTSTN(PTMX) .STRESS (PTMX) .DILAT(PTMX) , 

*  PRBSRV(PTMX) .SORRAD(PTMX) .SORPAR(PTMX) .SORVLP(PTMX) , 

*  IPDIST(PTMX) 

C 

IF  (IABORT. EQ.l)  RETURN 
C 

C*  load  master  arrays 
DO  20  I  =  l.NDIST 
DO  10  J  =  l.NPTS 
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S0RRAD((I-1)*NPTS+J)  =  RADIUS (I, J) 

S0RPAR( (I-1)*NPTS+J)  =  NUMPAR(I,J) 

SORVLP ( (I- 1)*NPTS+J)  =  VOLPAR(I.J) 

IPDIST((I-1)*NPTS+J)  =  I 
10  CONTINUE 
20  CONTINUE 
C 

C*  sort  master  arrays  in  ascending  order 
NTOT  =  NDIST*NPTS 

CALL  S0RT3 (NTOT , SORRAD , SORPAR , SORVLP , IPDIST) 

C 

C*  sort  master  arrays  in  descending  order 
DO  30  I  =  l,NDIST*NPTS/2 
ATMP  =  SORRAD (I) 

SORRAD (I)  =  SORRAD (NDIST*NPTS-I+1) 

SORRAD (NDIST*NPTS-I+i)  =  ATMP 
BTMP  =  SORPAR(I) 

SORPAR(I)  =  SORPAR (NDIST*NPTS-I+1) 

SORPAR (NDIST*NPTS-I+1)  =  BTMP 
CTMP  =  SORVLP (I) 

SORVLP (I)  =  SORVLP (NDIST*NPTS-I+1) 

SORVLP (NDIST*NPTS-I+1)  =  CTMP 
ITMP  =  IPDIST(I) 

IPDIST(I)  =  IPDIST (NDIST*NPTS-I+1) 

IPDIST (NDIST*NPTS-I+1)  =  ITMP 
30  CONTINUE 
C 

RETURN 

END 

C 

C 

SUBROUTINE  VOLFRC (NDIST , NPTS , VOLSMP , I ABORT) 

C====  calculates  dA/dc,  net  Vf,  net  Vv  and  probility  of  survival  for 
C  given  particle  radius.  Note:  net  Vf  is  based  on  total  sample  vol. 

C  Prob  of  surv  is  based  on  numbers  of  particles. 

C 

C  set  PTMX  =  NTDIS*GSMX 
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C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
INTEGER  GSMX , PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750.NTDIS  =  3) 

PARAMETER  (PI  =  3.1415927) 

COMMON  /DEBUG/  NUMPAR(NTDIS , GSMX) ,VOLPAR(NTDIS , GSMX) , 

*  NETVF (PTMX) , NETVV (PTMX) ,DADC(PTMX) , NPARTL (NTDIS) 

COMMON  /DIST/  RADAVG(NTDIS) .LOGSTD(NTDIS) .VLFRFO(NTDIS) , 

*  VLFRVO (NTDIS) 

COMMON  /RESULT/  CRTSTN(PTMX) , STRESS (PTMX) ,DILAT(PTMX) , 

*  PRBSRV(PTMX) .SORRAD(PTMX) .SORPAR(PTMX) .SORVLP(PTMX) , 

*  IPDIST(PTMX) 

C 

IF  (IABORT.EQ.l)  RETURN 
C 

C==  calculate  total  volume  fraction  filler  and  void 
NETVF ( 1 )  =  0 
NETVV (1)  =  0 
DO  10  I  =  1 ,NDIST 

NETVF (1)  =  NETVF ( 1 ) +VLFRF 0(1) 

NETVV (1)  =  NETVV ( 1 ) +VLFRV 0(1) 

10  CONTINUE 

PRBSRV(l)  =1.0 
C 

C==  calculate  net  Vf  and  Vv,  dA/dc  and  Prob  surv.  array  index  offset 

C  by  1  to  leave  room  for  initial  undebonded  state 

C 

SRVNUM  =  0 
C 

C==  find  total  number  of  particles 
DO  30  ICNT  =  l.NDIST 

TLNUMP  =  TLNUMP  +  NPARTL (ICNT) 

30  CONTINUE 
C 

DO  20  JCNT  =  2 ,NDIST*NPTS+1 

VLFTOT  =  VLFTOT  -  SORVLP( JCNT-1) 

NETVF (JCNT)  =  NETVF( JCNT-1)  -  SORVLP (JCNT-1) /VOLSMP 
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NETVV(JCNT)  =  NETVV(JCNT-l)  +  S0RVLP(JCNT-1)/V0LSMP 
SRVNUM  =  SRVNUM  +  SORPAR(JCNT-l) 

PRBSRV(JCNT)  =  (TLNUMP  -  SRVNUM) /TLNUMP 
DADC(JCNT)  =  -6. 0*V0LSMP/S0RRAD( JCNT-i) 

20  CONTINUE 
C 

RETURN 

END 

C 

C 

SUBROUTINE  MTPRP (CONCI , CONCV , ICNT , FDBND , YMULT , IKIND , IMORI , IPO IS , 

*  IABORT) 

C===s=  program  for  calculating  composite  modulus  based  on  Mori-Tanaka. 

C  FDBND=fraction  debond  for  orthotropic  properties  in  loading 
C  direction,  IKIND=use  inclusion  or  void  properties  in  calc  of 
C  Wv  matrix,  IMORI =type  of  particle  interaction  used  0=none, 

C  l=inclusion,  2=inclusion  and  void  or  vacuole, IPOIS=type  of 

C  debond  properties  O=orthotropic, l=isotropic 
C 

REAL  IDENT , K , KCMP , MAG 

PARAMETER  (GSMX=250 , PTMX=750 , NTDIS=3) 

COMMON  /PROPB/  Cll(PTMX) ,C12(PTMX) ,C21(PTMX) ,C22(PTMX) ,C23(PTMX) , 

*  ECMP(PTMX) ,POISC(PTMX) 

DIMENSION  CAVG(3 ,3) 

C 

IFCIABORT.EQ.1)  RETURN 
C 

IF ( ICNT. EQ.l) THEN 

CALL  CALCIO (IABORT) 

CALL  CALCCV (FDBND, IPOIS, IABORT) 

CALL  CMPRPO (IKIND, IMORI, IABORT) 

ELSE 
END  IF 
C 

CALL  CMPRP (CONCI , CONCV .YMULT , CAVG , IABORT) 

Cll (ICNT)=CAVG(1 , 1) 

C12(ICNT)=CAVG(i ,2) 
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C21(ICNT)=CAVG(2,1) 

C22(ICNT)=CAVG(2 ,2) 

C23(ICNT)=CAVG(2 ,3) 

ECMP(ICNT)=C11(ICNT)-2.0*C12(ICNT)*C21(ICNT)/(C22(ICNT)+C23(ICNT))  . 
P0ISC(ICNT)=C21(ICNT) / (C22( ICNT) +C23 (ICNT) ) 

C 

RETURN 

END 

C 

C 

SUBROUTINE  CRIT ( ICNT , VOLSMP , GAMM , PRESS , I ABORT) 

C«==  calculates  critical  strain  based  on  current  properties  because 
C  the  energy  balance  requires  the  input  work  to  equal  the  energy 
C  released  by  surface  creation  and  the  internal  energy  stored 

C  after  debonding  has  taken  place. 

C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
REAL  IDENT , K , KCMP , KTMP , MAG 
INTEGER  GSMX, PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750, NTDIS  =  3) 

PARAMETER  (TOL  =  IE-18) 

COMMON  /GAUS/  Z(GSMX) .RADIUS (NTDIS, GSMX) ,PROB(NTDIS ,GSMX) 

COMMON  /DEBUG/  NUMPAR (NTDIS , GSMX) ,VOLPAR(NTDIS ,GSMX) , 

*  NETVF (PTMX) .NETVV (PTMX) .DADC(PTMX) , NPARTL (NTDIS) 

COMMON  /PROPB/  Cil(PTMX) ,C12(PTMX) ,C21(PTMX) ,C22(PTMX) ,C23(PTMX) , 

*  ECMP (PTMX) ,P0ISC (PTMX) 

COMMON  /RESULT/  CRTSTN (PTMX) .STRESS (PTMX) ,DILAT (PTMX) , 

*  PRBSRV (PTMX) .SORRAD (PTMX) .SORPAR(PTMX) .SORVLP(PTMX) , 

*  IPDIST(PTMX) 

C 

IF  (IAB0RT.EQ.1)  RETURN 
C 

CONV  =  1.0E+3 
C 

DC  =  NETVF (ICNT) -NETVF ( ICNT- 1) 

IF(ABS(DC) .LT.TOL)  DC  =  -TOL 
C 
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TC12=C12(ICNT) 

TC21=C21(ICNT) 

TC22=C22(ICNT) 

TC23=C23(ICNT) 

DC11=(C11 (ICNT) -Cll (ICNT-1) ) /DC 
DC12=(C12(ICNT)-C12(ICNT-1))/DC 
DC21= (C21 (ICNT) -C21 (ICNT-1) ) /DC 
DC22=(C22(ICNT)-C22(ICNT-1))/DC 
DC23= (C23 (ICNT) -C23 (ICNT- 1) ) /DC 
C 

AQUAD  =  -DCil+2 .0* ( (TC22+TC23) * (TC2i*DC12+TC12*DC21)- 
*  (TC12*TC21* (DC22+DC23) ) ) / (TC22+TC23) **2 

CQUAD  =  C0NV*2*GAMM*DADC(ICNT) /VOLSMP 
C 

IF  ( AQUAD. GE.O)  THEN 

WRITE  (6, ’(A)’)  ’  SBR  CRIT:  square  root  term  is  negative.’ 
CRTSTN (ICNT)  =  CRTSTN( ICNT-1) 

ELSE 

CRTSTN (ICNT)  =  SQRT (CQUAD/AQUAD) 

END  IF 
C 

C==  for  debugging 

C  WRITE (6 , ’ ( A , 13) ’ )  ’ICNT:  ’ ,ICNT 

C  WRITE (6, ’ (A,E13.6,A,E13.6,A,E13.6) ’  ) 

C  *  ’  DE/DC:  ’ ,DCMPE, ’  GC:  ’.GAMM,’  DA/DC:  ’.DADC(ICNT) 

C  WRITE (6, ’ (A,E13.6,A,E13.6,A,E13.6) ’ ) 

C  *  ’  AQUAD:  ’.AQUAD,’  CQUAD:  ’.CQUAD,’  STN :  ’ .CRTSTN (ICNT-1) 

RETURN 
END 
C 
C 

SUBROUTINE  CAL VAL (ICNT, PRESS, DILATO, IABORT) 

C====  calculates  stress  and  dilatation  at  critical  strain 
C  properties  used  are  those  before  debonding  takes  place 
C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
REAL  IDENT , K , KCMP , KTMP , MAG 
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INTEGER  GSMX, PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750, NTDIS  =  3) 

COMMON  /PROPB/  Cll(PTMX) ,C12(PTMX) ,C21(PTMX) ,C22(PTMX) ,C23(PTMX) , 

*  ECMP(PTMX) ,POISC(PTMX) 

COMMON  /RESULT/  CRTSTN(PTMX) , STRESS (PTMX) ,DILAT(PTMX) , 

*  PRBSRV (PTMX) .SORRAD (PTMX) .SORPAR(PTMX) ,SORVLP(PTMX) , 

*  IPDIST(PTMX) 

C 

IF  (IABORT.EQ . 1)  RETURN 
C 

TC21  =  C21(ICNT-1) 

TC22  =  C22(ICNT-1) 

TC23  =  C23(ICNT-1) 

ETMP  =  ECMP(ICNT-l) 

IF  (ICNT.EQ .2)  DILATO  =  PRESS*0 
C 

STRESS ( I CNT)  =  ETMP*CRTSTN(ICNT) 

C==  change  stress  values  to  MPa 

STRESS (I CNT)  =  STRESS (ICNT)/1.0E6 

DILAT(ICNT)  =  ( 1- (2 . 0*TC2i/ (TC22+TC23) ) ) *CRTSTN(ICNT) -DILATO 
C 

RETURN 

END 

C 

C 

BLOCK  DATA  INIT 

C====  initialize  all  variables  and  arrays  used  in  program 
C  check  PTMX  if  NTDIS  or  GSMX  are  changed. 

C  PTMX  =  NTDIS*GSMX 

C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
REAL  IDENT , K , KCMP , MAG 
INTEGER  GSMX.PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750, NTDIS  =  3) 

COMMON  /GAUS/  Z(GSMX) .RADIUS (NTDIS, GSMX) ,PROB (NTDIS, GSMX) 

COMMON  /DEBUG/  NUMPAR(NTDIS, GSMX) ,VOLPAR(NTDIS, GSMX) , 

*  NETVF (PTMX) .NETVV(PTMX) .DADC(PTMX) .NPARTL (NTDIS) 
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COMMON  /DIST/  RADAVG(NTDIS) ,LOGSTD(NTDIS) ,VLFRFO(NTDIS) , 

*  VLFRVO(NTDIS) 

COMMON  /MATRA/  BETA(2) ,WI(3,3) ,WV(3,3) ,IDENT(3,3) 

COMMON  /MATRB/  S(3,3) ,CA(3,3) ,CB(3,3) ,CE(3,3) ,CF(3,3) 

COMMON  /PROPA/  K(3) ,G(3) ,E(3) ,P0IS(3) ,01(3,3) ,CV(3,3) ,C0(3, 3) 
COMMON  /PROPB/  Cil(PTMX) ,C12(PTMX) ,C21(PTMX) ,C22(PTMX) ,C23(PTMX) , 

*  ECMP(PTMX) .POISC(PTMX) 

COMMON  /RESULT/  CRTSTN(PTMX) .STRESS (PTMX) .DILAT(PTMX) , 

*  PRBSRV (PTMX) ,SORRAD(PTMX) .SORPAR(PTMX) .SORVLP(PTMX) , 

*  IPDIST(PTMX) 

C 

DATA  Z  /GSMX*0/  RADIUS  /PTMX*0/  PROB  /PTMX*0/ 

DATA  NUMPAR  /PTMX*0/  VOLPAR  /PTMX*0/  NETVF  /PTMX*0/  NETVV  /PTMX*0/ 

*  DADC  /PTMX*0/  NPARTL  /NTDIS*0/ 

DATA  RADAVG  /NTDIS*0/  LOGSTD  /NTDIS*0/  VLFRFO  /NTDIS*0/  VLFRVO  / 

*  NTDIS*0/ 

DATA  BETA  /2*0/  HI  /9*0/  WV  /9*0/  IDENT  / 1,0, 0,0, 1,0, 0,0,1/ 

DATA  S  /9*0/  CA  /9*0/  CB  /9*0/  CE  /9*0/  CF  /9*0/ 

DATA  K  /3*0/  G  /3*0 /  E  /3*0/  POIS  /3*0/  Cl  /9*0/  CV  /9*0/  CO  /9*0/ 
DATA  Cll  /PTMX*0/  C12  /PTMX*0/  C21  /PTMX*0/  C22  /PTMX*0/ 

*  C23  /PTMX*0/  ECMP  /PTMX*0/  POISC  /PTMX*0/ 

DATA  CRTSTN  /PTMX*0/  STRESS  /PTMX*0/  DILAT  /PTMX*0/ 

*  PRBSRV  /PTMX+O/  SORRAD  /PTMX*0/  SORPAR  /PTMX*0/  SORVLP  /PTMX*0/ 

*  IPDIST  /PTMX*0/ 

C 

END 

C 

C 

SUBROUTINE  QSIMP (FUNC , A ,B , S) 

C==*=  used  for  integration  of  gaussian  curve  in  sbr  GAUSS,  obtained 
C  from  Numerical  Recipes,  W.H.  Press,  Cambridge,  1988. 

C 

EXTERNAL  FUNC 

PARAMETER  (EPS  =  1.E-6.JMAX  =  20) 

OST  =  -1.E30 
OS  =  -1.E30 
DO  10  J  =  1 , JMAX 


P501263.PDF  [Page:  97  of  122] 


UNCLASSIFIED 

B.19 


CALL  TRAPZD(FUNC,A,B,ST,J) 

S  ■  (4.*ST-0ST)/3 . 

IF  (ABS(S-OS) .LT.EPS*ABS(OS))  RETURN 
OS  =  S 
OST  =  ST 
10  CONTINUE 

PAUSE  ’Too  many  steps:  SBR  QSIMP’ 

END 

C 

C 

SUBROUTINE  TRAPZD(FUNC,A,B,S,N) 

C====  used  for  integration  of  gaussian  curve  in  sbr  QSIMP  which  is 
C  called  from  sbr  GAUSS,  obtained  from  Numerical  Recipes,  W.H. 

C  Press,  Cambridge,  1988. 

C 

EXTERNAL  FUNC 
IF  (N.EQ.l)  THEN 

S  =  0 .5*(B-A)*(FUNC(A)+FUNC(B) ) 

IT  =  1 
ELSE 

TNM  =  IT 
DEL  =  (B-A)/TNM 
X  =  A+0 . 5*DEL 
SUM  =  0. 

DO  10  J  =  1 ,IT 

SUM  =  SUM+FUNC(X) 

X  =  X+DEL 
10  CONTINUE 

S  =  0.5*(S+(B-A)*SUM/TNM) 

IT  -  2*IT 
END  IF 
RETURN 
END 
C 
C 

SUBROUTINE  S0RT3(N,RA,RB,RC,IRD) 

C====  sorting  routine  from  Numerical  Recipes,  W.H.  Press,  Cambridge, 
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C  1988.  sorts  in  ascending  order  array  RA  and  moves  elements  in 
C  arrays  RB,  RC  and  IRD  at  the  same  time. 

C 

DIMENSION  RA(N) ,RB(N) ,RC(N) ,IRD(N) 

L  =  N/2+1 
IR  =  N 
10  CONTINUE 

IF  (L.GT.l)  THEN 
L  =  L-l 
RRA  =  RACL) 

RRB  =  RB(L) 

RRC  =  RC(L) 

IRRD  =  IRD(L) 

ELSE 

RRA  =  RA(IR) 

RRB  =  RB(IR) 

RRC  =  RC(IR) 

IRRD  =  IRD(IR) 

RA(IR)  =  RA(1) 

RB(IR)  =  RBCl) 

RC(IR)  =  RC(1) 

IRD(IR)  =  IRD(l) 

IR  =  IR-1 
IF  (IR.EQ.l)  THEN 
RA(1)  =  RRA 
RB(1)  =  RRB 
RC(1)  =  RRC 
IRD Cl)  =  IRRD 
RETURN 
ENDIF 
END  IF 
I  =  L 
J  =  L+L 

20  IF  (J.LE.IR)  THEN 

IF  (J.LT.IR)  THEN 

IF  (RA(J) .LT.RACJ+l))  J  =  J+l 


ENDIF 
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IF  (RRA.LT.RA(J))  THEN 
RA(I)  =  RA(J) 

RB(I)  =  RB(J) 

RC(I)  =  RC(J) 

IRD(I)  =  IRD(J) 

I  =  J 
J  =  J+J 
ELSE 

J  =  IR+1 
ENDIF 
GOTO  20 
ENDIF 

RA(I)  =  RRA 
RB(I)  =  RRB 
RC(I)  =  RRC 
IRDCI)  =  IRRD 
GOTO  10 
END 
C 
C 

SUBROUTINE  CALCIO(IABORT) 

C====  calculate  the  property  matrix  for  inclusion  and  matrix, 

C  isotropic  relations 
C 

REAL  IDENT , K , KCMP , MAG 

COMMON  /PROPA/  K(3) ,G(3) ,E(3) ,P0IS(3) ,CI(3,3) ,CV(3,3) ,C0(3,3) 
C 

IF (IABORT .EQ . 1)  RETURN 
C 

K(l)  =  (2.O*G(l)*(l+P0Isa)))/(3.O*(i.O-2.O*P0IS(i))) 
E(1)=G(1)*(2 .0*(1+P0IS(1) )) 

K(3)=(2 . 0*G(3) *(i+P0IS(3) ) ) / (3 .0* (1 .0-2 .0*POIS(3) ) ) 
E(3)=G(3)*(2.0*(1+POIS(3))) 

Ci=K(l)+(4.0/3.0)*G(l) 

C2=K(l)-(2.0/3.0)*G(l) 

C3=K(3)+(4.0/3.0)*G(3) 

C4=K(3)-(2.0/3.0)*G(3) 
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DO  11  1=1,3 
DO  21  J=1 ,3 
CI(I , J)=C2 
C0(I, J)=C4 

IF(I.EQ.J)  CI(I,J)=C1 
IF(I.EQ.J)  C0(I, J)=C3 
21  CONTINUE 
11  CONTINUE 
C 

RETURN 

END 

C 

C 

SUBROUTINE  CALCCV (FDBND , IPOIS , I ABORT) 

C====  calculate  the  property  matrix  for  debonded  particle, 

C  orthotropic  relations,  FDBND  is  debond  fraction  for  vacuole 
C  IPOIS  determines  whether  orthotropic  or  isotropic 
C 

REAL  IDENT , K , KCMP , HAG 

COMMON  /PROPA/  K(3) ,G(3) ,E(3) ,P0IS(3) ,CI(3,3) ,CV(3,3) ,C0(3,3) 
C 

IF(IABORT.EQ.i)  RETURN 
C 

IF(K(2).NE.O.AND.G(2).NE.O)  THEN 

P0ISC2)  =  (3.0*K(2)-2.0*G(2))/ (2.0*(3.0*K(2)+G(2))) 

EC2)  =  9 . 0*K(2) *G(2) / (3 . 0*K(2)+G(2) ) 

ELSE 

E(2)  =  0.0 
P0IS(2)  =  0.0 

END  IF 
C 

PC0N=REAL (IPOIS) 

DETM  =  1-P0IS(2)**2-PC0N*2*(P0IS(2)**2+P0IS(2)**3) 

CV(i , 1)=(FDBND*E(2)*(1-P0IS(2)**2))/DETM 

CV(1 ,2)=(FDBND*E(2)*(P0IS(2)+P0IS(2)**2))/DETM 

CV(1 ,3)=CV(1 ,2) 

CV(2, 1)=(E(2)*PCON*POIS(2)*(1+POIS(2)))/DETM 
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CV (2 , 2)= (E(2) *(l-PC0N*P0IS(2)**2) ) /DETM 
CV(2 ,3)-(E(2)*(P0IS(2)+PC0N*P0IS(2) **2) ) /DETM 
CV(3, 1)=CV(2, 1) 

CV(3,2)=CV(2,3) 

CV(3,3)=CV(2,2) 

C 

RETURN 

END 

G 

C 

SUBROUTINE  CMPRPO (IKIND , IMORI , IABORT) 

C====  calculate  constants  in  composite  equation 
C 

REAL  IDENT , K , KCMP , MAG 

COMMON  /MATRB/  S(3,3) ,CA(3,3) ,CB(3,3) ,CE(3,3) ,CF(3,3) 

COMMON  /PROPA/  K(3) ,G(3) ,E(3),P0IS(3) ,CI(3,3),CV(3,3) ,C0(3,3) 
DIMENSION  CTEMPA (3,3), CTEMPB (3,3) 

C 

IF ( IABORT. EQ.l)  RETURN 
C 

CALL  CALCW( IKIND, IMORI, I ABORT) 

CALL  CALCS (IABORT) 

CALL  SUB (CTEMPA, Cl, CO) 

CALL  INVERT (CTEMPB , CTEMPA , IABORT) 

CALL  MULT (CA, CTEMPB, CO) 

CALL  SUB (CTEMPA, CV, CO) 

CALL  INVERT (CTEMPB , CTEMPA , IABORT) 

CALL  MULT (CB, CTEMPB, CO) 

CALL  ADD(CE,S,CB) 

CALL  ADD(CF,S,CA) 

C 

RETURN 

END 

C 

SUBROUTINE  CMPRP (CONCI , CONCV , YMULT , CAVG , IABORT) 

C==  calculate  composite  properties,  ITYPE  identifies  inclusion 
C  or  void 
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REAL  IDENT , K , KCMP , MAG 

COMMON  /MATRA/  BETA(2) ,WI(3,3) ,WV(3,3) ,IDENT(3,3) 

COMMON  /MATRB/  S(3,3) ,CA(3 ,3) ,CB(3,3) ,CE(3,3) ,CF(3,3) 

COMMON  /PROPA/  K(3) ,G(3) ,E(3) ,P0IS(3) ,CI(3 ,3) ,CV(3,3) ,C0(3,3) 
DIMENSION  CC(3,3) ,CD(3,3) ,CG(3,3) ,CH(3,3) 

DIMENS ION  CTEMPA (3,3), CTEMPB (3,3), CTEMPC ( 3 , 3 ) , C AVG (3 , 3 ) 

C 

C==  calculate  phase-dependent  components  of  composite  equation 

C  calculate  first  half 

C  calculate  phase-i  components 

I TYPE  =  1 

CALL  GAMMA ( CG , CONCI , I TYPE , YMULT , I ABORT) 

CALL  SUB (CTEMPA, IDENT, S) 

CALL  SUB (CTEMPB, CTEMPA, CG) 

DO  11  1=1,3 
DO  21  J=1 ,3 

CC (I , J) =CONCI*CTEMPB (I , J) 

21  CONTINUE 

11  CONTINUE 

C  calculate  phase-v  components 

ITYPE  =  2 

CALL  GAMMA (CH , CONCV , ITYPE , YMULT , I ABORT) 

CALL  SUB (CTEMPA, IDENT, S) 

CALL  SUB ( CTEMPB , CTEMPA , CH ) 

DO  31  1=1,3 
DO  41  J=1 ,3 

CD (I , J)=CONCV*CTEMPB(I , J) 

4:1  CONTINUE 

31  CONTINUE 

CALL  INVERT (CTEMPA, CE.IABORT) 

CALL  MULT ( CTEMPB , CD , CTEMPA ) 

CALL  MULT (CTEMPA, CTEMPB, CF) 

C  combine  phase-i  and  phase-v  components 
CALL  ADD (CTEMPB, CTEMPA, CA) 

CALL  ADD (CTEMPA, CTEMPB, S) 

CALL  ADD (CTEMPB, CTEMPA, CC) 

CALL  INVERT ( CTEMPA , CTEMPB , I ABORT) 
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CALL  MULT (CTEMPB ,CG, CTEMPA) 

DO  51  1=1,3 
DO  61  J=1 ,3 

CTEMPC (I , J)=CONCI*CTEMPB(I , J) 

61  CONTINUE 
51  CONTINUE 

C  calculate  second  half 

CALL  INVERT (CTEMPA ,CF , IABORT) 

CALL  MULT ( CTEMPB , CC , CTEMPA ) 

CALL  MULT (CTEMPA, CTEMPB, CE) 

C  combine  phase-i  and  phase-v  components 

CALL  ADD (CTEMPB, CTEMPA, CB) 

CALL  ADD (CTEMPA, CTEMPB, S) 

CALL  ADD (CTEMPB , CTEMPA , CD ) 

CALL  INVERT (CTEMPA .CTEMPB , IABORT) 

CALL  MULT ( CTEMPB, CH, CTEMPA) 

DO  71  1=1,3 
DO  81  J=1 ,3 

CTEMPA (I , J)=CONCV* CTEMPB (I, J) 

81  CONTINUE 
71  CONTINUE 

C=  combine  all  components 

CALL  ADD (CTEMPB, CTEMPC, CTEMPA) 

CALL  ADD (CTEMPA, CTEMPB, IDENT) 

CALL  MULT (CAVG, CO, CTEMPA) 

C 

RETURN 

END 

C 

C 

SUBROUTINE  CALCW (IKIND , IMORI , IABORT) 

C====  calculate  correction  matrices  WI  and  WV  and  BETA  for 
C  use  in  sbr  GAMMA,  IKIND  determines  inclusion  or  void  for  vacuole 
C  IMORI  determines  if  correction  matrix  used,  O=none, i=inclusion 
C  2=inclusion  and  void 
C 


REAL  IDENT, K.KCMP, MAG 
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REAL  KTEMP , KMAT 

COMMON  /MATRA/  BETA(2) ,WI(3,3) ,WV(3,3) ,IDENT(3 ,3) 

COMMON  /MATRB/  S(3,3) ,CA(3,3) ,CB(3,3) ,CE(3,3) ,CF(3,3) 

COMMON  /PROPA/  K(3) ,G(3) ,E(3) ,P0IS(3) ,CI(3,3) ,CV(3,3) ,C0(3,3) 
C 

IF(IABORT.Eq.l)  RETURN 
C 

POISM  *  PQIS(3) 

KMAT  =  K(3) 

GMAT  =  G(3) 

C 

KTEMP=K(1) 

GTEMP=G (1) 

DO  11  INCL=1 ,2 

IF (INCL . EQ . 2 . AND . IKIND . EQ . 0)THEN 
GTEMP=0 . 0 
KTEMP=0 . 0 
ELSE 

GTEMP=G(INCL) 

KTEMP=K(INCL) 

END  IF 
C 

ALPHA  =  2.0*(5 .0*POISM-1)+10 .O*(i-POISM)* 

*  (KMAT/ (KTEMP-KMAT) -GMAT/ (GTEMP-GMAT) ) 

BETA (INCL)  =  2 . 0* (4 . 0-5 . 0*POISM) +15 . 0* ( 1-POISM) * 

*  (GMAT/ (GTEMP-GMAT)) 

ZETA1  =  12 . 0* ( 13 . O+POISM-14 . 0*P0ISM**2) - (96 . 0* ALPHA/ 

*  (3 . 0*ALPHA+2 . 0*BETA (INCL)  ) ) * ( 1-2 . 0*P0ISM)  *  ( 1+POISM) 
ZETA2  =  6.0*(25.0-34.0*POISM+22.0*POISM**2)-(36.0*ALPHA/ 

*  (3 . 0*ALPHA+2 . 0*BETA (INCL) ) ) * ( 1-2 . 0*P0ISM) * ( 1+POISM) 

C 

DO  21  1=1,3 
DO  31  J=1 ,3 

IF (INCL . EQ . 1 . AND . IMORI . NE . 0) THEN 
WI(I, J)=ZETAi 

IF(I.EQ.J)  WI(I,J)=ZETA1+2*ZETA2 
ELSEIF (INCL . EQ . 1 . AND . IMORI . EQ . 0) THEN 
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WI(I,J)=0.0 

ELSEIF (INCL . EQ . 2 . AND . IMORI . EQ . 2) THEN 

WV(I,J)=ZETA1 

IF(I.EQ.J)  WV(I,J)=ZETA1+2*ZETA2 
ELSE 

WV(I,J)=0.0 
END  IF 

31  CONTINUE 

21  CONTINUE 
11  CONTINUE 
C 

RETURN 

END 

C 

C 

SUBROUTINE  CALCS (IABORT) 

C====  calculate  Eshelby  matrices  SI  and  SV 

COMMON  /MATRB/  S(3,3) ,CA(3,3) ,CB(3,3) ,CE(3,3) ,CF(3,3) 

COMMON  /PROPA/  K(3) ,G(3) ,E(3) ,P0IS(3) ,CI(3,3) ,CV(3,3) ,C0(3,3) 
C 

IF ( IABORT. EQ.i)  RETURN 
C 

POISM  =  P0ISC3) 

SDET  =  15.0*(1-PDISM) 

C 

51  =  5 ,0*POISM-1 

52  =  4.0-5.0*POISM 
C 

DO  11  1=1,3 
DO  21  J=1 ,3 

S(I , J)=S1/SDET 

IF(I.EQ.J)  SCI, J)=(S1+2.0*S2) /SDET 
21  CONTINUE 
11  CONTINUE 
C 

RETURN 

END 
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C 

SUBROUTINE  GAMMA (A ,C0NC , ITYPE ,YMULT , I ABORT) 

C====  calculate  correction  matrix  A  given  inclusion  I  and  its 
C  concentration  CONC,  Y  depends  on  microstructural  features 
C  ITYPE  identifies  inclusion  or  void 
REAL  IDENT , K , KCMP , MAG 

COMMON  /MATRA/  BETA(2) ,WI(3,3) ,WV(3,3) ,IDENT(3,3) 
DIMENSION  A (3, 3) 

C 

IF (IAB0RT.EQ . 1)  RETURN 
C 

Y=YMULT* ( i-CONC) /24 . 0 

MAG  =  5.0*CONC*Y/ (4.0*BETA(ITYPE)**2) 

DO  11  1=1,3 
DO  21  J-1,3 

IF (ITYPE . EQ . 1)  A (I , J)=IDENT(I , J)+MAG*WI(I , J) 
IFCITYPE.EQ.2)  A(I , J)=IDENT(I , J)+MAG*WV(I, J) 

21  CONTINUE 
11  CONTINUE 
C 

RETURN 

END 

C 

C 

SUBROUTINE  ADD(C,A,B) 

C====  subroutine  for  adding  two  square  matrices  C=A+B 
DIMENSION  A(3,3),B(3,3),C(3,3) 

C 

DO  11  1=1,3 
DO  21  J=1 ,3 

C(I,J)  =  A(I , J)+B(I , J) 

21  CONTINUE 
11  CONTINUE 
C 

RETURN 

END 
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C 

C 

SUBROUTINE  SUB(C,A,B) 

C====  subroutine  for  adding  two  square  matrices  C=A-B 
DIMENSION  A(3,3) ,B(3,3) ,C(3,3) 

C 

DO  11  1=1,3 
DO  21  J-1,3 

C(I,J)  =  A(I,J)-B(I,J) 

21  CONTINUE 
11  CONTINUE 
C 

RETURN 

END 

C 

C 

SUBROUTINE  MULT(C,A,B) 

C====  subroutine  for  multiplying  two  square  matrices  C=A.B 
DIMENSION  AC3.3) ,B(3,3) ,C(3,3) 

C 

DO  11  1=1,3 
DO  21  J=1 ,3 
C(I,J)  =  0 
DO  31  K=1 ,3 

C(I,J)  =  C(I,J)+A(I,K)*B(K,J) 

31  CONTINUE 

21  CONTINUE 
11  CONTINUE 
C 

RETURN 

END 

C 

C 

SUBROUTINE  INVERT (AI, A, I ABORT) 

C====  subroutine  used  for  inverting  matrix  A  to  give  AI 
DIMENSION  A(3,3) ,AI(3,3) 


C 
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IF(IABORT.EQ.i)  RETURN 
C 

DETA  =  -(A(1,3)*A(2,2)*A(3,1))  +  A(1 ,2)*A(2,3) *A(3, 1) 

*  +  A(1,3)*A(2,1)*A(3,2)  -  A(1,1)*A(2,3)*A(3,2) 

*  -  A(1,2)*A(2,1)*A(3,3)  +  A(1,1)*A(2,2)*A(3,3) 

C 

IF (DETA. NE.O) THEN 

AI(l,i)  =  (-(A(2,3)*A(3,2))  +  A(2 ,2)*A(3 ,3) ) /DETA 
AI(1, 2)  =  (A(1,3)*A(3,2)  -  A(1,2)*A(3,3))/DETA 
AI(1, 3)  =  (-(A(1,3)*A(2,2))  +  A(1 ,2)*A(2,3) )/DETA 
AI(2,i)  =  (A(2,3)*A(3,1)  -  A(2,1)*A(3,3))/DETA 
AIC2.2)  =  (-(A(1,3)*A(3,1))  +  A(1,1)*A(3,3))/DETA 
AI(2,3)  =  (A(l ,3)*A(2, 1)  -  A(1,1)*A(2,3))/DETA 
AIC3.1)  =  (-(A(2,2)*A(3,1))  +  A(2,1)*A(3,2))/DETA 
AI(3,2)  =  (A(i ,2) *A(3, 1)  -  A ( 1 , 1)*A (3,2) ) /DETA 
AI(3,3)  =  (-(A(l ,2)*A(2, 1) )  +  A(l,i)*A(2,2))/DETA 
ELSE 

I ABORT  =  1 

WRITE(6, ’ (A) ’)  ’SBR  INVERT:  indeterminant  matrix’ 

END  IF 
C 

RETURN 

END 

C 

C 

SUBROUTINE  GAUWRT (NDIST , NPTS , I ABORT) 

C====  write  out  cumulative  distribution  data. 

C  for  some  reason,  cannot  print  out  PROBs  correctly  using 
C  F  format,  numbers  end  up  getting  multiplied  by  ten. 

C 

C  set  PTMX  =  NTDIS*GSMX 
C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
INTEGER  GSMX.PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750,NTDIS  =  3) 

COMMON  /GAUS/  Z (GSMX) .RADIUS (NTDIS, GSMX) ,PROB(NTDIS, GSMX) 
C 
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IF  (IAB0RT.EQ . 1)  RETURN 

WRITE  (6, ’(A)’)  ’  Writing  GAUSS . DAT ’ 

C 

OPEN  (UNIT=7 , FILE= ’ _GAUSS . DAT ’ ,F0RM=’ FORMATTED’ ,STATUS= ’UNKNOWN’ ) 
WRITE  (7,5000) 

DO  10  IPTS  =  1 , NPTS 

WRITE  (7,5100)  Z(IPTS), (PROB(IDIST, IPTS) ,RADIUS(IDIST, IPTS), 

*  IDIST  =  l.NDIST) 

10  CONTINUE 

CLOSE  (7) 

C 

RETURN 

5000  FORMAT  (’  Z  Pr  Radius (mm)  Pr 

*  Radius (mm)  Pr  Radius (mm)’) 

5100  FORMAT  (1X,F6.3,6(3X,0PE13.6)) 

END 

C 

C 

SUBROUTINE  HSTWRT (NDIST , NPTS , IABORT) 

C====  write  out  histogram  and  tracking  data 
C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
INTEGER  GSMX,PTMX 

PARAMETER  (GSMX  =  250,PTMX  =  750.NTDIS  =  3) 

COMMON  /GAUS/  Z (GSMX) .RADIUS (NTDIS, GSMX) ,PROB(NTDIS, GSMX) 

COMMON  /DEBUG/  NUMPAR(NTDIS , GSMX) ,VOLPAR(NTDIS , GSMX) , 

*  NETVF (PTMX) .NETVV(PTMX) .DADC(PTMX) .NPARTL (NTDIS) 

C 

IF  (I ABORT. EQ . 1)  RETURN 

WRITE  (6, ’(A)’)  ’  Writing  HISTO .DAT’ 

C 

OPEN  (UNIT=7,FILE=’ _HISTO.DAT’ ,FORM=’ FORMATTED’ , STATUS=’ UNKNOWN ’ ) 
NPRTOT  =  0 
VOLTOT  =  0 

DO  30  IDIST  =  l.NDIST 
DO  40  IHST  =  l.NPTS 

NPRTOT  =  NPRTOT  +  NUMPAR(IDIST.IHST) 
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V0LT0T  =  V0LT0T  +  VOLPAR(IDIST,IHST) 

40  CONTINUE 
30  CONTINUE 
C 

CUMVOL  =0.0 
DO  20  IDIST  =  l.NDIST 
WRITE  (7,5000) 

DO  10  IHST  =  l.NPTS 

PERNPR  =  100*REAL(NUMPAR(IDIST, IHST)) /REAL (NPRTOT) 

PERVOL  =  100*VOLPAR(IDIST , IHST) /VOLTOT 

CUMVOL  =  CUMVOL  +  PERVOL 

WRITE  (7,5100)  IHST, RADIUS (IDIST, IHST), 

*  ALOG10(NUMPAR(IDIST, IHST)), VOLPAR(IDIST, IHST) , PERNPR, 

*  PERVOL, PROB (IDIST, IHST) , CUMVOL 
10  CONTINUE 

20  CONTINUE 
CLOSE  (7) 

C 

RETURN 

5000  FORMAT  ( ’  Point  avg  R(mm)  log  #  part.  volume  (mm3)  '/.no. 

♦particles  '/,  part. volume  cum.  prob.  cum.  vol.’) 

5100  FORMAT  (2X,I3,3X,7(1PE13.6,2X)) 

END 

C 

C 

SUBROUTINE  STRWRT (NDIST , NPTS , VOLSMP , GAMM , FDBND , YMULT , IKIND , IMORI , 

*  IPOIS , PRESS , DILATO , FILNM , I ABORT) 

C====  write  out  stress  and  dilatation  results  versus  critical  strain 
C  include  probability  survival,  radius,  no.  particles  and 

C  distribution  info. 

C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
REAL  IDENT , K , KCMP , MAG 
INTEGER  GSMX , PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750.NTDIS  =  3) 

COMMON  /DIST/  RADAVG(NTDIS) , LOGSTD (NTDIS) ,VLFRFO(NTDIS) , 

*  VLFRVO (NTDIS) 
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COMMON  /DEBUG/  NUMPAR(NTDIS ,GSMX) ,VOLPAR(NTDIS ,GSMX) , 

*  NETVF(PTMX) .NETVV(PTMX) ,DADC(PTMX) .NPARTL(NTDIS) 

COMMON  /PROPA/  K(3) ,G(3) ,E(3) ,P0IS(3) , Cl (3,3) ,CV(3,3) ,C0(3,3) 
COMMON  /PROPB/  Cll(PTMX) ,C12(PTMX) ,C2l(PTMX) ,C22(PTMX) ,C23(PTMX) , 

*  ECMP(PTMX) .POISC(PTMX) 

COMMON  /RESULT/  CRTSTN(PTMX) , STRESS (PTMX) .DILAT(PTMX) , 

*  PRBSRV (PTMX) ,SORRAD(PTMX) .SORPAR(PTMX) .SORVLP(PTMX) , 

*  IPDIST(PTMX) 

CHARACTER  FILNM*8 

C 

IF  (I ABORT. EQ . 1)  RETURN 
C 

IF  (FILNM . EQ . 5  DEFAULT  * )  FILNM  =  ’ .STRESS’ 

WRITE  (6, ’ C/,A,A8,A) ’)  5  Writing  to  5 , FILNM, » .DAT’ 

C 

OPEN  (UNIT=7 , FILE=FILNM// » . DAT ’ , STATUS= 5  UNKNOWN 5 ) 

WRITE  (7,5000) 

DO  10  I  *  l.NDIST 

WRITE  (7 ,  ’ (1X,I1 ,4(3X ,0PE11 .4) )  ’ )  I ,RADAVG(I) ,LOGSTD(I)  , 

*  VLFRFO(I) ,VLFRV0(I) 
iO  CONTINUE 

C 

WRITE  (7,5100)  G(3) ,G(1) ,P0IS(3) ,P0IS(1) ,G(2) ,K(2) 

WRITE  (7,5200)  VOLSMP,FDBND,YMULT,IKIND,IMORI,IPOIS 
WRITE  (7,5300)  PRESS, GAMM,DILATO 
WRITE  (7,5400) 

DO  20  I  =  1 ,NDIST*NPTS+1 
ETMP  =  ECMP(I)/1E6 

WRITE  (7,5600)  I ,CRTSTN(I) .STRESS (I) ,DILAT(I) , PRBSRV (I) ,ETMP , 

*  POISC(I)  ,NETVF (I)  ,NETW(I)  .IPDIST(I-l) 

20  CONTINUE 

C 

CLOSE  (7) 

RETURN 

5000  FORMAT  (»  #  avg  Rad(um)  std  dev  Vf  Vv’) 

5100  FORMAT  (’  Gm(Pa)=» ,0PE11 .4, »  Gf (Pa)=’ .0PE11.4, 5  vm=’ ,0PE11 .4, 

*  »  vf=’ ,0PE11.4,’  Gv (Pa) = ’ , 0PE1 1.4,’  Kv=’ ,0PE11 .4) 
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5200  FORMAT  (»  V(mm3)  *»,0PE11.4,»  frac  dbnd=’ ,0PEii.4, »  Y-mult=’, 

*  0PE11 .4, 5  w-type=’ ,13, »  m-type=’ ,13,  ’  v-type=’,I3) 

5300  FORMAT  (»  P0(Pa)=’ .0PE11.4, »  Gc(Pa-m)=’ ,0PE11 .4, ’  (dV/V)0=» ,0PEii 

*  .4) 

5400  FORMAT  (’  Point  crit  strn  stress (MPa)  dV/V  Prlsurv 

*  E_c(MPa)  Poisson  V_f  V_v  diet’) 

5600  FORMAT  (1X,I3,3X,8(1PE11.4,2X),1X,I1) 

C 

END 

C 

G 

SUBROUTINE  DBGWRT (NDIST, NPTS , I ABORT) 

C===;=  write  out  additional  data  for  debugging  purposes 
C 

C  set  PTMX  =  NTDIS*GSMX 
C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
REAL  IDENT , K , KCMP , MAG 
INTEGER  GSMX.PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750,NTDIS  =  3) 

COMMON  /GAUS/  Z (GSMX) .RADIUS (NTDIS, GSMX) ,PROB(NTDIS, GSMX) 

COMMON  /DEBUG/  NUMPAR(NTDIS, GSMX) ,VOLPAR(NTDIS, GSMX) , 

*  NETVF (PTMX) , NETVV (PTMX) .DADC(PTMX) .NPARTL (NTDIS) 

COMMON  /PROPB/  Cl 1 (PTMX) , C12 (PTMX) ,C21 (PTMX) ,C22(PTMX) ,C23(PTMX) , 

*  ECMP (PTMX), POISC (PTMX) 

COMMON  /RESULT/  CRTSTN(PTMX) .STRESS (PTMX) .DILAT(PTMX) , 

*  PRBSRV(PTMX) .SORRAD(PTMX) .SORPAR(PTMX) .SORVLP(PTMX) , 

*  IPDIST(PTMX) 

C 

IF  (IABORT.EQ . 1)  RETURN 

WRITE  (6, ’ (A) 5 )  ’  Writing  DEBUG.DAT’ 

C 

OPEN  (UNIT=7 , FILE= ’ _DEBUG . DAT  * , FORM= ’ FORMATTED ’ , STATUS= 5  UNKNOWN ’ ) 
WRITE  (7,5000) 

DO  10  IHST  =  1 ,NDIST*NPTS+i 

WRITE  (7,5100)  IHST, NETVF (IHST) .NETVV (IHST) .DADC(IHST) , 

*  Cll(IHST) ,C12(IHST) ,C21(IHST) ,C22(IHST) ,C23(IHST) , 
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*  PRBSRV (IHST) 

10  CONTINUE 

CLOSE  (7) 

RETURN 

5000  FORMAT  ( ’  Point  net  Vf 

net  Vv 

dA/dc 

*  Cll  C12 

*  Pr I surv ’ ) 

C21 

C22 

C23 

5100  FORMAT  (2X,I3,3X,9(0PE13.6,2X)) 

END 

C 

G 

SUBROUTINE  DBGRAT (NDIST , NPTS , IABORT) 

C====  write  out  additional  data  for  debugging  purposes 
C  along  with  stress-strain  data  outputs  radius  and  the  factor 
C  SQRT(RAD*dG/dc)  to  look  at  its  relationship  with  crit.  strain 
C 

C  set  PTMX  =  NTDIS*GSMX 
C 

REAL  LOGSTD , NPARTL , NUMPAR ,NETVF , NETVV 
REAL  IDENT , K , KCMP , KTMP , MAG 
INTEGER  GSMX , PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750.NTDIS  =  3) 

COMMON  /GAUS/  Z(GSMX) , RAD IUS(NTDIS, GSMX) .PROB(NTDIS.GSMX) 

COMMON  /DEBUG/  NUMPAR(NTDIS , GSMX) ,VOLPAR(NTDIS , GSMX) , 

*  NETVF (PTMX) .NETVV(PTMX) ,DADC(PTMX) .NPARTL (NTDIS) 

COMMON  /PROPB/  Cll(PTMX) ,C12(PTMX) ,C21(PTMX) ,C22(PTMX) ,C23(PTMX) , 

*  ECMP(PTMX) .POISC(PTMX) 

COMMON  /RESULT/  CRTSTN (PTMX) .STRESS (PTMX) .DILAT (PTMX) , 

*  PRBSRV (PTMX) .SORRAD (PTMX) .SORPAR(PTMX) .SORVLP(PTMX) , 

*  IPDIST(PTMX) 

C 

IF  (IABORT.EQ . 1)  RETURN 

WRITE  (6,’ (A)’)  »  Writing  DERAT . DAT ’ 

C 

OPEN  (UNIT=7 , FILE= 5  _DERAT . DAT ’ , FORM= 5  FORMATTED ’ , STATUS* ’ UNKNOWN ’ ) 
WRITE  (7.5000) 
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DO  10  IHST  =  1 ,NDIST*NPTS+i 
ETMP  =  ECMP(IHST) 

EAD  =  0.0 
DNETF  =0.0 
DNETV  =0.0 
DETMP  =0.0 

IF (IHST. GT. 1)  RAD  =  SORRAD(IHST-l) 

IF(IHST.GT.l)  DNETF  =  ABS(NETVF(IHST)  -  NETVF(IHST-l) ) 
IF(IHST.GT.l)  DNETV  =  ABS(NETVV(IHST)  -  NETVV(IHST-l) ) 

IF (IHST. GT. 1)  DETMP  =  ABS(ECMP(IHST)  -  ECMP(IHST-l)) 

FACT  =  SQRT (RAD*DETMP) 

WRITE  (7,5100)  IHST, CRTSTN (IHST) .STRESS (IHST) , RAD, PRBSRV(IHST) , 

*  ETMP, POISC (IHST) .DNETF, DNETV, FACT 
10  CONTINUE 

CLOSE  (7) 

C 

RETURN 

5000  FORMAT  (’  Point  crit  strn  stress (MPa)  Avg  r(mm)  Prlsurv 

*  E_c(MPa)  Poisson  dV_f  dV_v  fact’) 

5100  FORMAT  (IX, I3,3X,9(1PE11 .4,2X) ) 

END 

C 

C 

INCLUDE  ’MSGRAPH .FOR’ 

C 

SUBROUTINE  CRVPLT(NDIST , NPTS , IABORT) 

C====  driver  routine  for  plotting  curve  on  screen,  keep  the 
C  INCLUDE  ’MSGRAPH. FOR’  with  this  module. 

C 

C  set  PTMX  =  NTDIS*GSMX 
C 

REAL  LOGSTD , NPARTL , NUMPAR , NETVF , NETVV 
INTEGER  GSMX , PTMX 

PARAMETER  (GSMX  =  250, PTMX  =  750.NTDIS  =  3) 

COMMON  /RESULT/  CRTSTN (PTMX) .STRESS (PTMX) .DILAT (PTMX) , 

*  PRBSRV (PTMX) .SORRAD (PTMX) .SORPAR(PTMX) .SORVLP(PTMX) , 

*  IPDIST(PTMX) 
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DIMENSION  X(PTMX) .Yl(PTMX) ,Y2(PTMX) 

CHARACTER  ANS*1 
C 

IF  (IAB0RT.EQ.1)  RETURN 
C 

10  CONTINUE 

WRITE  (6,’(/,A)’)  ’  Graph  results  on  screen?  (Y/N) ’ 
READ  (5 , 5 (Al) ’ )  ANS 
C 

NTOT  =  NDIST*NPTS+1 
DO  20  ICNT  =  l.NTOT 

X(ICNT)  =  CRTSTN(ICNT) 

Y1 (ICNT)  =  STRESS (ICNT) 

Y2(ICNT)  =  DILAT(ICNT) 

20  CONTINUE 
C 

IF  (ANS.EQ.’Y’)  THEN 

WRITE  (6, ’(A)’)  ’  Strain,  Stress  and  dV/V  end  pts’ 

READ  (5,*)  XEND , YSEND , YDEND 

CALL  GRAF (NTOT , X , Y1 , Y2 , XEND , YSEND , YDEND) 

ELSE 
END  IF 
C 

IF  (ANS.EO. ’Y’)  GOTO  10 
C 

RETURN 

END 

C 

C 

SUBROUTINE  GRAF(N,X,Y1 ,Y2 , XEND .YSEND .YDEND) 

C 

C  set  PTMX  =  NTDIS*GSMX 
C 

INTEGER  PTMX 

PARAMETER  (GSMX  =  250. PTMX  =  750.NTDIS  =  3) 

DIMENSION  X(PTMX) .Yl(PTMX) ,Y2(PTMX) 

C 
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INCLUDE  ’ GRFDEF .FOR’ 

C 

CALL  VIDEO (MAXX , MAXY , NOGRAF ) 

IF  (NOGRAF. EQ.O)  THEN 
C 

CALL  VWPORT (MAXX, MAXY) 

C 

XBEG  =  0 
YBEG  =  0 

CALL  WINDOW (XBEG , YBEG , XEND , YSEND) 

ICURV  =  1 
XLAB  =  ’strain’ 

YLAB  =  ’strs  (MPa)’ 

CALL  ATTRIB (ICURV, ILNCOL,ILNSTY) 

CALL  LABELS (ICURV , ILNCOL , XL AB , YLAB , XBEG , YBEG , XEND , YSEND) 
C 

ICURV  =  1 

CALL  ATTRIB (ICURV, ILNCOL, ILNSTY) 

CALL  CURVE(X,Y1,N, ILNCOL, ILNSTY) 

CALL  WINDOW (XBEG , YBEG , XEND , YDEND) 

ICURV  =  3 
XLAB  =  ’strain’ 

YLAB  =  ’dV/V’ 

CALL  ATTRIB (ICURV, ILNCOL, ILNSTY) 

CALL  LABELS (ICURV , ILNCOL , XLAB , YLAB , XBEG , YBEG , XEND , YDEND) 
ICURV  =  3 

CALL  ATTRIB (ICURV, ILNCOL, ILNSTY) 

CALL  CURVE(X,Y2,N, ILNCOL, ILNSTY) 

C 

CALL  ENDGRFO 
ELSE 

WRITE  (6, ’(A)’)  ’  SBR  GRAF:  problem  with  graphics’ 

ENDIF 

C 

RETURN 

END 
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